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Chapter 1 

introduction 



In physics the "system" is the small part of the universe on which we have some control. 
In other words, it is the part of "interest" that we measure and/or manipulate. Knowing 
the system Hamiltonian, TCs, we can calculate its time evolution and find the future state 
of the system, as well as telling its past (how it reached the current state). The system of 
interest can be, for instance, an electron, an atom, a piece of metal, a cup of coffee, or the 
solar system. 

Everyday experience teaches us that any system interacts with its environment. For 
example, we all know that a coffee thermos in a cold winter morning will eventually 
cool down, no matter how sophisticated the flask is. Nobody is astonished by this: we 
cannot fully isolate the coffee from the outside. Admittedly, in many cases the effects of 
the environment can be neglected, and the dynamics is then determined by the system 
Hamiltonian alone. In other cases, however, we must take into account that our system is 
not isolated, but coupled to its surroundings [HE]. 

Either because we do not know the evolution of the environment (on which we have 
not control), or because we only look at the system dynamics, we do not have complete 
info on the whole (system plus bath). Think about classical mechanics, for example; one 
point in the system's phase space may correspond to several configurations of the total 
system. If we cannot fully determine the current state, we are unable to reconstruct the 
past or predict the future evolution from TLs alone. On the other hand, the total phase 
space of system plus environment is very large. As a result the Poincare recurrence time 
(to return to the initial state) becomes infinite for most practical purposes: we never find 
the coffee warming up after some time (irreversibility). 

In this context the formalism of open systems is important in that we need not change 
the current theoretical frameworks to account for irreversibility. In addition, this formal- 
ism provides some clues for the difficulty in finding quantum superpositions in everyday 
macroscopic life. The theory saves face. 

In the following pages we study the effects of the environment on the dynamics of the 
system of interest. Although we start surveying the classical limit, our work will focus on 
the quantum domain (non-relativistic). A constant throughout the manuscript will be to 
relate the results obtained with their classical limits (appropriately defined). 
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Although the formahsm of the first chapters is general, the theory will be applied 
later on to systems with one degree of freedom. Examples from solid-state physics where 
the calculations can be applied are the phase difference across a Josephson junction, or 
collective excitations in lattices, like solitons or fluxons pll3]. In these systems openness 
is important, while they can be approximately described by a few collective degrees of 
freedom. We will also tackle spin systems, like quantum superparamagnets or ensembles 
of 1/2-spins (two-level systems) [HE]- 

As for the spirit of this manuscript, we have chosen not to repeat here the details 
of the calculations included in the articles we have written in the last three years. The 
reader interested in those details can consult our publications, which are cited as PAPER I, 
II,. . .In any case the thesis is quite self-contained, and can be followed without recourse 
to the articles (at least that was our aim). 

The text is organized as follows. We begin reviewing the classical limit in Chap. [2l 
We introduce the bath-of-oscillators formalism for the environment, obtaining reduced dy- 
namical equations (Langevin & Fokker-Planck equations). Then we quantize and proceed 
to derive dynamical equations for the reduced density matrix (quantum master equations; 
Chap. [3]). Chapter m is devoted to the thermal equilibrium properties of open quantum 
systems. Then in Chapter Owe introduce the phase-space formalism and apply it to quan- 
tum master equations, giving a natural link to the classical equations. This first part is 
general, and we close it by introducing in Chap. El the tools we will employ later on in 
explicit calculations: linear-response theory & the continued-fraction technique to solve 
master equations. The last two chapters (the longest) are devoted to the application of the 
general formalism to two problems of interest. In chapter [7] we study the Bloch electron 
(a particle in a periodic potential) coupled to a dissipative bath. Finally, in Chap. [8] we 
study the equilibrium and dynamical properties of quantum superparamagnets. 
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Chapter 2 



classical theory: 

irreversibility and quantization 



In this chapter we place the problem in a historical context, beginning with the 19th cen- 
tury discovery of Brownian motion. Then we introduce the bath-of-oscillators formalism 
and address in this frame the problems of irreversibility and quantization. We conclude 
introducing the specific systems we will study in this thesis. 

2.1 history of Brownian motion [6, Chap. 1] 

In 1827 Father Brown observes under his microscope the irregular motions executed by 
pollen grains when suspended in a fluid (figure [27T|) . Such a motion disagreed with physi- 
cal laws known at the time, where particles followed smooth and predictable trajectories 
according to Newton laws (we left aside chaos, unknown then). Due to his background, 
he was a botanist. Brown thought that the animated motion observed might be a mani- 
festation of life. Soon he would be discarding such an explanation. 




Figure 2.1: Possible trajectory observed by Brown. 
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It wasn't until 1905 when Einstein, in his annus mirabilis, explains in a satisfactory 
way the phenomenon, rightly called Brownian motion^ Einstein's approach is based on 
two main points: (i) the irregular motion is caused by the impact of the fluid's molecules 
on the colloidal grain (ii) the motion of those molecules is so intricate that the grain motion 
can only be described statistically. From these assumptions an equation for W{x, t) could 
be derived, giving the probability of the pollen grain being around x at time t (a sort of 
Liouville equation, but in configuration space). 

dW{x,t) ^^ d'^W^x.t) ^ 
dt dx'^ 

This is a diffusion equation where D is related with the standard deviation of the grain 
position {x{tf) - (x(0)2) = 2Dt. 

A few years later Langevin proposed an explanation, a model, which using his own 
words was "infinitely more simple" than that of Einstein & Smoluchowsky. He assumed 
that the force due to the fluid can be decomposed into: (i') a systematic part (viscous 
force) and (ii') a fluctuating force, accounting for the unpredictability of the collisions. 
Incorporating these forces into Newton's equation, Langevin wrote: 

d?x dx ^, , , , 

= -"^7^ + e(i) (2.2) 

Here 7 is the viscous damping parameter and ^{t) the fluctuating force, on which he only 
assumed that it can take random values both positive and negative. This was the first 
example of a stochastic differential equation [HJ [7] . 

Except for the short time regime (inertial), the treatments of Einstein and Langevin 
(more complete) lead to equivalent results. Both equation (12. ip and ()2.2p describe the 
pollen grain dynamics. The Liouville and Newton time evolutions, respectively, get mod- 
ified due to the interaction with the fluid, in such a way that the effective equations for 
the grain display irreversible dynamics. 



2.2 problems with irreversibility and quantization 

At this point we could raise the following question. Assuming that, though certainly 
complex, both the pollen grain and the fluid molecules obey Newton's reversible laws, 
how comes one can arrive at irreversible equations? This is the irreversibility paradox [2]: 

• Loschmidt objection: equations (j2.ip and (j2.2p are not invariant under time-reversal. 

To this, a second complain was added: El 

^ In 1906 Smoluchowsky 's paper is published, putting forward the same explanation of Brownian motion. 
We would like to dispel an urban legend attributing the name to the irregular and unpredictable motions 
on stage of the late James Brown. 

^ These objections were first raised against Boltzmann's kinetic equation [S]. 
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• Zermelo's paradox: If we start from a configuration witli many pollen grains at the 
center of the fluid, we know that after some time they will spread throughout the 
volume, according to the diffusion equation ()2.ip . That is, we will never find that 
at some later time all grains return to the center, the initial configuration. This is 
at odds with Poincare's theorem, which asserts that the trajectories of a bounded 
system (in phase space) will pass arbitrarily close to the initial state, after a time 
called the recurrence time. 

Finally, along with the irreversibility /recurrence issues, we find the problem of quan- 
tization. How could we quantize an equation that is irreversible, or stochastic? In general 
phenomenological quantization schemes (say, the counterpart of plugging —7 x in Newton 
equation) can lead to violations of basic issues (like the very normalization of the state, 
or commutation relations) H With the formalism introduced in next section we will try to 
answer both questions. 

2.3 bath of oscillators formalism 

Here we discuss the rigorization of the diffusion and Langevin equations. This will allow 
us to address the issues of irreversibility and quantization. 

As physicists we feel comfortable when we see a Hamiltonian. Therefore we formally 
write: 



This total Hamiltonian is made of the system Hamiltonian Tig, the environment T^b (the 
bath), and the interaction between them Tish- Naturally, we will have to give content 
to each part of Wtot- This looks like a though task though. On top of the problem of 
modelling the system itself (think of nuclear physics), now we have to model everything in 
its surroundings. Fortunately, in many cases of interest one has the following simplifying 
conditions: 

(i) The environment "feels" only weakly the presence of the system. Loosely speaking, 
the bath is big, a macroscopic object (but not necessarily classical); the liquid carrier 
in Brown's set up. The fluid's properties change very little when a pollen grain, or 
a handful, are suspended on it. Naturally, the converse does not hold. The grain is 
indeed affected by the bath. 

(ii) The environment dynamics is faster that the characteristic time scales of the system. 
For example, due to their microscopic nature the fluid molecules move around so fast 
that they look like random to us. 

Under these conditions the bath can be modelized as a set of harmonic oscillators 
linearly coupled to the system |TT] : 



^ Several attempts of direct quantization of dissipative equations are critically discussed in [9] and in 



(2.3) 





[3 Sec. 2.2]. 
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Here {x,p) are the system canonical variables, not necessarily coordinates and velocities. 
Note that we can extract the interaction as ~ QaF{x,p). Thus, the coupling is linear in 
the bath coordinates Qa-, whereas it can be arbitrary in the system variables F{x,p). 

When extracting the coupling ~ QaF{x,p), one also gets a term oc F'^{x,p) which 
only depends on the system variables. This counter-term cancels the renormalization 
of 7is produced by a purely linear ~ QaF (generally speaking, an arbitrary coupling 
induces fluctuations, dissipation and renormalizations). So what we are doing is assuming 
that Tis is the Hamiltonian to which we have access experimentally, where all relevant 
renormalizations have already been included (see [TO] and what follows). 

A further case where Hamiltonian (12.41) can be used is: 



(iii) The bath oscillators provides an "exact" description of the environment. Think 
about the phonons in a solid, or the normal modes of the electromagnetic field, the 
photons, or less evident cases like electron-hole excitations. 

Well, once we have introduced our working Hamiltonian, we would like to obtain 
dynamical equations for the (sub)system degrees of freedom. To this end we consider a 
dynamical variable A{x^p) depending only on the system, and write Hamilton's equations 
for it: dA/dt = {A,7itot} with { , } the Poisson bracket and the total Hamiltonian ([27 
Writing similarly the equations for the bath variables, Qa and Pa, we have |12j 



a a 

dQa D dPa 2 r\ (c^ n\ 

= Pa, = -UJaQa - UaF . (2.7) 

The first line includes the contribution of the system dynamics in the absence of interaction 
(first term), the part from the renormalization terms (middle), and the coupling to the 
bath coordinates (last). The second line is simply the equations of motion for the bath 
oscillators (as forced by the system). 

Here we start to appreciate the advantages of model (j2.4p . as we know how to solve 
Eqs. (12. 7p . They are just scores of non-interacting and forced oscillators, with F{x,p) 
playing the role of a driving. This was made possible by the coupling being linear in 
the bath coordinates. Now, from undergraduate physics we know that the solution of a 
harmonic oscillator forced by F{t) reads: 

Qa = Ql- — f dr s\n[LOa{t - t)]F{t) . (2.8) 

Jto 



* The most general Hamiltonian linear in the bath would be 

Htot=H,{x,p) + (Pc + ^;«Gc(x,p))%lj^(q„ + ^Fc,{x,p)y (2.5) 

where one accounts for a possible coupling with the oscillator momenta, and for different F and G for each 
mode. However, to illustrate the formalism is enough with p.4|l . The general (|2.5p brings little new, while 
it makes the calculations cumbersome. 
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where Q\ = Qa{tQ) cos[u)a{t — to)] + Pa{to)/^aS^'^[^a{t — to)] is the general solution of 
the homogeneous equation (free oscillators, without F), while the last term is a particular 
solution of the inhomogeneous equation. 

Plugging now the solution (j2.8p in the system equation (j2.6p . and integrating it by 
parts, one finds {j3 = l/k-QT) 



A A 

— = {A,n.} + {A,F} 



f{t)+(5 j\TC{t-T)'^{T) (2. 



,2 



We see that the renormalization term in p. 60 was cancelled by the term — -^^{A, F^} 
arising from the integration by parts. This was indeed the reason behind including ~ F'^ 
in the starting Hamiltonian ()2.4p . It ensures that the renormalization (if physical) was 
already included in Ji^^ and not counted twice; then "Hg is the Hamiltonian (the levels) to 
which we have access in experiments. As for the other two terms in (12. 9p . we are going to 
see they describe fluctuations and relaxation with 

/(t) = ^t.„Q^; C(t-T) =fcBr^^cosK(t-r)] (2.10) 

a a ^ 

We have extracted /? from C{t — t) so that it matches the correlator of /(t) [Eq. (|2.1ip 
below], in correspondence with the quantum definition (chapter [3]) . 

From the equations above we begin to grasp the structure of a Langevin equation like 
()2.2p . We would have the "random force" /(t), which is just the free evolution of the 
bath oscillators (see Fig. 12. 2p . while the memory integral in Eq. (j2.9p would correspond 
to the viscous damping. The later originates from the inhomogeneous solution in (12. Sp . 
and brings back the dynamics at previous times (a sort of back-reaction on the system of 
its previous action on the bath). As for the fluctuations, recall that we do not have full 
control on the bath, whereas /(t) includes the initial conditions (5o(to), -Po(to) for each 
oscillator. It seems reasonable to draw them from an equilibrium state at the initial timeU 
Being harmonic oscillators, the distribution of positions and momenta would be Gaussian, 
so we just need to specify the first two moments of /(t) 

{fit)) = (2.11) 
(/(t)/(r)) = C{t-r). (2.12) 

Here we see the link mentioned between the integral kernel and the force correlator. The 
latter can be written in terms of the bath spectral density: 

J{u;) = 7Ty^^6{io-iOa) (2.13) 

a 

in the following way [13] 

C(r) = 2fcBT r^:Mcosu;r. (2.14) 
Jo 

Now we are ready to give a couple of examples to give content and fix ideas. 



® This does not mean that the bath is at equihbrium for t > to, since the system induces dynamics on 
the bath through the interaction; see the Qa evolution 
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Figure 2.2: f{t) obtained from ()2.10p summing over 1000 oscillators with initial conditions 
drawn from a Gaussian distribution. The force really looks random. 



2.3.1 examples 

Brownian particle 



-X I.e., 



Let us consider the Hamiltonian of mechanics TCs = p'^ /2m + V{x), and F{x,p) 
the coupling is linear in the system coordinate as well). Assume now that the correlator 
C{t) is a function strongly peaked around r = 0, expressing that the system loses memory 
quickly (Markovian limit). Quantitatively C(t) = 2mk-QT^ 5{t), corresponding to J{uj) = 
m^uj, with 7 a measure of the coupling strengthlf] Setting first A = x and then A = p i 
the equation of motion ()2.9p one gets 



m 



dx 
dt 



dm 



dp 
dt 



dm 



+ fit) 



dx 



(2.15) 



dp dt dx 

This equation reduces to the one proposed by Langevin when Tig = p^ /2m (that is, for 
the free particle). 

A stochastic equation with Gaussian f{t) and ^-correlation {f{t)f{T)) oc 6{t — t), 
can be converted into an equation for the probability distribution W (Fokker-Planck 
equation; see e.g. [21 [7]). For the mechanical problem ()2.15p one obtains the Klein-Kramers 
equation 



dtW{x,p,t) = - p/mdx + V'dp + -/dp{p + mkBTdp) W{x,p,t) 



(2.16) 



This Fokker-Planck equation for a particle in a potential is a generalization of the diffusion 
equation (f2T]l . [l| 



If a = 1 the bath is called 



Ohmic (Markovian case) and the cases a ^ 1 are called super/sub-Ohmic. The choice a = 1 if often done 



The spectral density (|2.14p is customarily assumed of the form J{u!) 

for simplicity; it is hard to solve non-local (integro-differential) equations! But if required by the physics 
of the problem, one cannot escape working with a 7^ 1, as we will do later on. 

^ Actually, Einstein Eq. (|2.ip follows when neglecting inertial effects, the so called overdamped regime. 
And it was for a free particle; in a potential the overdamped equation for the distribution W{x,t) = 
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classical Brownian spin/dipole 



Fortunately the equation of motion (12. 9|) is valid for any Poisson bracket, so we can 
handle problems more general than the Brownian particle. Given that in this thesis 
we will study simple spin Hamiltonians, we are going to discuss first the classical case 
TCs = 7is{mx,rny,mz) with the angular- momentum Poisson algebra 



m 



(sin cos </>, sin 9 sin cos 6) ; 



{mi,mj} = geijkrrik 



(2.18) 



Here 6 and (p are the polar and azimuthal angles on the sphere. The constant g is the gy- 
romagnetic ratio related with the spin number by hS = fJ-s/g, with ^ub Bohr's magnetron. 
In this case we set A = rrii in (|2.9|) and we get (Ohmic bath) H 



1 dm 

5^ 



in X 



Bef + m 



1 

s 



m X A( m X B 



cf 



B 



cf 



dm 



(2.20) 



This equation has the structure of the phenomenological equation proposed by Landau- 
Lifshitz in 1935. Here {d/dm)i = d/drrii is the gradient in m-space, while the noise and 
the damping include: 



m = fit) 



OF 
dm 



A 



jk 



A 



OF dF 
drrij druk 



(2.211 



We have used an Ohmic J{u>) = Xlj, with A the coupling strength; the analogue of the 
damping 7 in particle problems (footnote [6]) . We will use two diff'erent letters for particles 
and spins following the standard notation in the literature (besides, 7 has dimensions of 
frequency 1/t whereas A is dimensionless) . 

As we did for the Brownian particle, we can also write a Fokker-Planck equation 
associated to the spin Langevin equation (j2.20p 



I dW 



d 
dm 



[m X Bef) — ^m X A 



m X ( B 



kBT 



d 
dm 



W 



(2.22) 



where now (d/dm)- is the divergence operator in m-space. Let us illustrate this with 
a simple case: coupling i*" ~ ^ linear in the spin variables plus axially symmetric 
7is{0), as in Tig oc —K cos^ 6 — BzCosO. Under these conditions S^Bef = and we can 
restrict ourselves to axially symmetric solutions d^W = 0. Then, introducing z = cos 9 
and 2rN = (3/{gX), equation (|2.22|) reduces to [T3]: 



2rN- 



dW 



dt 



d_ 

dz 



{l-z') 



d ^dn, 

dz dz 



W 



J dpW{x,p;t) is called the Smoluchowsky equation [14] 

'ydtWix,t) = {kBTdl,-d^V)W{x,t) 



(2.23) 



(2.17) 



Note the useful relation: 



{mi,/(m)} = -g m x 



dfjm) 
dm 



(2.19) 



valid for any function of the spin variables rrii. It follows from the Poisson bracket (|2.18p and the definition 
of the vector product x . 
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This is nothing but Debye's equation for dipolar relaxation (possibly the first Fokker- 
Planck equation describing rotational Brownian motion). 

2.4 revisiting the paradoxes of irreversibility 

Now we are in a position permitting to address the two objections of section 12.21 

• First, let us tackle the time-reversal objection due to Loschmidt. It is clear that the 
Hamilton equations (|2.6p and (|2.7|) are reversible. However, the generalized Langevin 
equation (12. 9p is completely equivalent to them. What happens is that one thinks 
of reversing the velocities of the system only. Therefore, we have a "practical" 
irreversibility: although the whole system is reversible, we cannot change the sign 
of the velocities of all degrees of freedom, including those of the bath (on which in 
practice we have little control). 

• Second, Zermelo's paradox. The fact that both the effective equations (and the 
experiments) show that the pollen grains are never brought back together, is just a 
consequence of the total phase space being very large. Again, we cannot overlook 
the bath. Indeed one can estimate the recurrence time and find that it diverges as 
we let the number of bath oscillators tend to infinity [U [16]. It is only for finite 
"baths" that such time remains finite. But in practical cases the recurrence time 
can be several times the estimated age of the universe. For practical purposes we 
can be sure that the grains will never return to the initial configuration. 

2.5 quantization 

We have just seen how to handle the irreversibility problem (get along with it, would 
be more correct). Let us move on to the issue of quantizing a system with Langevin or 
Fokker-Planck dynamics. The idea is to consider the Hamiltonian (|2.4|) as the problem to 
be quantized. Then we would just need to transform the variables {x,p) & {Qa,Pa) into 
Hilbertian operators and proceed by the book. 

Actually, the main derivation of Sec. 12.31 can be followed analogously in the quantum 
case; "simply" replace Poisson brackets { , } by commutators [ , ], with A an operator 
in the Heisenberg picture. This is how one derives quantum Langevin equations \17] . 
However, these are not as useful/generic as their classical counterparts; now both x and 
the force f{t) are operators, and calculations to-the-end exist only in simple cases \18\ I19j 
On top of this, there are no simple recipes (as those in the classical case), to convert a 
Langevin equation into an equation for the distribution (Fokker-Planck like). 

To bypass these problems one usually resorts to quantum master equations. Recall 
that the bath is not under control, so it is out of question any preparation of pure states 
(of system plus bath). Then one needs the density matrix formulation (more general) 
to compute averages {A), our main interest. If A acts only on the Hilbert space of our 
system, we can introduce the reduced density matrix g. Then using the properties of 
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Figure 2.3: Left: energy bands for the cosine potential (j2.26p with K = 200 (the quantum- 
ness parameter ~ Sq/H). Right: ratchet potential ()2.27p with r=0.44 and bands plotted 
for K = (see the text). 



partial tracing, we get {A) just as a trace over the system: 

{A) = Trs+b(^£'tot) = Trs(A£)) ; q := Trb(£>tot) 



(2.24) 



Here s + b indicates the original trace over the whole system+bath using the total density 
matrix ^ptot- 

Much as in the classical case the distribution W{x,p) involved the system variables 
only, now we have g playing this role. In this sense the equations of motion for g will be 
the analogue of the Fokker-Planck equations. Those equations (next chapter) generalize 
the Von Neumann equation for the closed evolution: 



dg 

'dt ~ 

This simply parallels the fact that the Hamiltonian part of the Fokker-Planck was the 
Liouville equation 

dtW = {ns,w} . 



2.6 systems studied in this thesis 

We are going to deal with both particle problems and spin Hamiltonians (translational and 
rotational problems, respectively). In the final applications we will consider Hamiltonians 
with one degree of freedom (but which could describe effectively composite systems like a 
Josephson junction or a molecular magnet). 

2.6.1 particle systems 

We will consider the standard Hamiltonian 



(2.25) 
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with [x,p] = ih. In particular we will address the problem of a quantum particle running 
in a periodic potential V{x) = V{x + L) plus a constant force (this is related with the 
Wannier-Stark problem [20]). We will study the dynamics in a cosine potential (Fig. I2.3] l: 

y(x) = Vbcosx (2.26) 

and in ratchet potentials (lacking inversion symmetry) 

V{x) = -Vo [smx + (r/2) sin(2x)] . (2.27) 

In figure [231 we have plotted the energy bands for these potentials (see PAPER I, Sec. 
7.3). We have introduced a dimensionless parameter 

K = 27ry ; (2.28) 

with Sq some relevant/characteristic action; for instance So = Eq/loq with £^0 the barrier 
height and loq the oscillation frequency at the bottom of the wells. The parameter K tells 
us how classical (or quantum) the system is. Letting — > 00 gives the classical limit; 
then the bands tend to a continuum (for K = 200 the levels of Fig. 12.31 are already very 
close) . 



2.6.2 spin systems 

In parallel with mechanical problems, we will consider simple spin Hamiltonians (Fig. 12. 4|) 

Hs = -DSl - B,S, . (2.29) 

Here [Si, S^] = isijkSk, with S'^\m) = S{S + l)\m) and Sz\m) = \m). The above is an effec- 
tive spin Hamiltonian [21j . The term —B^Sz is the Zeeman coupling to the external field, 
while D is the anisotropy constant, of spin-orbit origin. When Z) = 0, Hamiltonian (j2.29p 
reduces to the often studied isotropic spin. The cases D ^ correspond to anisotropic 
spins, with "easy axes" and "easy plane" anisotropy, respectively. Later, in chapter [8] we 
will focus on D > in (|2.29p . which is the minimal model for superparamagnets [22]. 

The classical limit will be taken letting S = {nB/g)/fi tend to infinity. Note the 
analogy between K for periodic potentials and S here. When taking the limit, it will 
be convenient to "normalize" the operators as Hs = —DS'^{Sz/ S)'^ — BzS{Sz/ S). Then 
we let S grow keeping constant DS^ and BzS, the amount of anisotropy and Zeeman 
energies. Mathematically, more and more levels are introduced, towards a continuum, 
keeping constant the barriers' height. 
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Figure 2.4: Energy levels of the spin Hamiltonian (|2.29|) for 5 = 5 and different anisotropy 
and field parameters. Top left: isotropic spin {D = 0). Top right: spin with D < at 
Bz = \D\; the spectrum has a single well structure. Bottom: double minima structure, 
D > 0, in the absence of applied field (left) and at the first crossing field = D (right). 
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Chapter 3 



dynamics of open quantum 
systems 



In this chapter we address the dynamics of open quantum systems. We first discuss 
the notion of quantum decoherence, through a solvable example. To handle more complex 
systems one can resort to quantum master equations, the quantum analogue of the Fokker- 
Planck equations of the previous chapter. Here we introduce those equations and give 
some examples. We discuss the approximations involved in their derivation, as well as 
other invoked in solving master equations. 

3.1 quantum decoherence 

In the previous chapter we discussed how irreversibility could emerge from reversible equa- 
tions. Now we would like to see how the classical world we are familiar with arises out of 
the microscopic quantum world it is made of. A proposed explanation for such emergence 
resorts to the coupling with the environment, in a irreversible process called decoherence 
(see |23) for an introduction to the trendy topic). It can be seen as a consequence of 
having the dynamics of very many degrees of freedom but only looking at some reduced 
dynamics. The notion can be introduced with a simple example, due to Zurek [24j, which 
is fully solvable. 

The question one aims to answer is: if everything around us, including ourselves, 
is made of atoms obeying the laws of quantum mechanics, how comes we do not find 
superpositions of states in the macroscopic world? The most famous/shocking of those 
superpositions is due to Schrodinger: 



Or others equally exotic like |^') oc | here) + | there), i.e., the possibility of having delocalized 
objects. To address these points let's go through the promised simple examplej^ 

^ It is worth recalling that the superposition (jSTTJ is diflterent from a statistical mixture (as Schrodinger 
himself noted). In the density-matrix language, a statistical mixture and a superposition are represented 




(3.1) 
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As usual, we will consider that the system interacts with many degrees of freedom, so 
that the total quantum Hamiltonian is Titot = + ^sb + ^b- For simplicity, we assume 
that both system and bath have discrete and non-degenerate spectra. In the basis of 
eigenstates of Tls and TCy, one can write 

T-Cs = '^6i\si){si\ ; Hh = ^ea\ea){ea\ ■ (3.3) 

i a 

Here le^) are the eigenstates of the bath oscillators, and \sj) those of the system. In this 
basis, the TCs^ in the bath-of-oscillators model (j2.4p can be written as 



'^7ia\si){si\ |eQ,)(eo| + ^ |Sj)(Sj/| \e-a){ea'\ 

ia ii'aa' 

To proceed we assume that (Jaiaa' ^ 7iai meaning that decoherence times are much 
shorter than relaxation processes. Then 

Wsb = ^^lia\\Si){Si\ \ea){ea\] , (3.4) 



which allows to solve analytically the problem. If an'aa' ^ 7ia we would have to resort to 
a quantum master equation (next section), and solve it numerically, most likely. 

Let us assume now that our system is in a superposition like (13. ip and let us study the 
dynamical evolution. Since we have to specify the state of the bath too, we assume for 
simplicity that it is a pure state (though everything can be repeated for the realistic case 
of a mixed bath). Then, the initial state reads: 

\'^{t = 0)) = Y,aj\sj) (^^ba\ea) (3.5) 

j a 

Now with the "Hgb considered, the time evolution at time t results in the entangled stateH 

l^'(t)) = ^aj6„e-i*(^^+^"+^^-)|sj) ® \ea) (3.6) 

Calculating now the reduced density matrix (I2.24p . that is, ^ = Trb(|\l'(t))(\l'(t)|) one finds 
in the basis 

Sjj{t) = (3-7) 
Qijit) = aia*e~-^'(-'^-'^hij{t); Zi,-(t) = ^ |6„| V^""* (3.8) 



by 



„,i,ed / 1/2 \ * / 1/2 1/2 , 

^ = ( 1/2 j ' ^ = U/2 1/2 ' • ^^-^^ 



The state (|3.5p is "separable", that is, it can be expressed as the tensor product of a state of the sys- 
tem's Hilbert space times a state in the bath's Hilbert space. On the contrary, (|3.6p is a hnear combination 
of separable states, and it is no longer possible to assign state vectors to the system and bath separately. 
These states are called entangled. 
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Figure 3.1: Several realizations of a random walk in 2 dimensions: z = |6a|e^'^''*. The 
total sum is constrained as X] |&aP = 1 (see text). 

Here just the coherent rotation due to the evolution under Tig, while 

Zij{t) is due to the interaction with the environment [uja = {^ia — lja)/K\. 

In general ufi 7^ and hence Zij{t) results to be the sum of many oscillating factors 
with different frequencies. Then at long enough times {zij{t)) 0, so that e~"^"* can be 
seen as a random phase. Therefore e~"^° becomes a random vector in the complex 
plane. The problem can thus be mapped onto the random walk in 2D, where each move 
is l^al^e"^^"* and ^ l^aP = 1 is kept fixed (see Fig. 13. ip . Then Zij{t) is a random variable 
with a probability distribution given by P{z) oc e~^^^ , where is the number of bath 
oscillators. Therefore at sufficiently long times, and considering the bath to be macroscopic 
(A^ — > oo), one writes 

z,j{t)^0. (3.9) 

Or in other words, the off-diagonal elements drop to zero and the original superposition 
X] j I ) is "destroyed" . In the sense that the reduced density matrix is indistinguishable 
from a statistical mixure with I diagonal elements. 

The time td, from which we can take Zij{t) = 0, is called the decoherence time. The 
loose of coherence is a consequence of looking only at the system. In fact, everything 
follows from the unitary evolution of the total quantum system; that is, the coherence is 
still out there, but now is diffused/shared between a large number of degrees of freedom. 
In a sense, this decoherence follows from the "non-locality" of quantum mechanics, as the 
state (13. 6p is fully entangled between system and bath. 

Finally, is worth to recall that every density matrix is always diagonalizable (as it is 
Hermitian) . Or conversely, given a diagonal density matrix, it can be seen as non-diagonal 
by some other basis. Then, the concept of decoherence seems to be "basis dependent". 
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What happens, decoherentists say, is that the interaction with the bath "determines" in 
which basis q will be diagonal; the so called pointer basis [23\ . It seems that Nature chose 
us to be localized objects, as well as being alive or death. To give some numbers, the 
decoherence time for a superposition of dust grains with a diameter 10~^ cm, set 1 cm 
apart is estimated to be around 10~^^ seconds [25] . That is, for most practical purposes, 
and for most macroscopic objects, quantum superpositions cannot be observedo 



3.2 derivation of master equations 

It will not be always possible to express Q{t) analytically, as in the previous section. Our 
goal here is to derive differential equations for g. In contrast to the classical case [Eq. (12. 9p ]. 
we do not know how to obtain in general an exact equation for the model (j2.4p . Therefore, 
we will be happy with an approximate equation. The main assumption involved will be 
to consider that the system-bath interaction is weak. 

The most usual techniques in deriving master equations are those of projection opera- 
tors [21 [28], path integral P [291 [301 [3l] , or cumulant expansions [3 [32]. For our purposes, 
it will suffice with a simpler derivation, based on plain quantum-mechanical perturbation 
theory. The approximations invoked will be clearly seen as we proceed. 

We will start working with gtot we would trace later on to obtain the effective dynamics 
for Q [Eq. (I2.24p ]. The evolution of gtot is given by 

Qtotit) = U{t; to)gtot{to)U+{t; to) ; Uit; to) = e-i«-*(*-*o) (3.10) 

where we have assumed that Htot is not time dependent. This is not essential though; 
we could consider that Wtot depends on t, and then the evolution operator would be a 
time-ordered exponential U{t;to) = T[exp{—i/h J dt'Wtot (*')}] ■ 

In what follows it will be convenient to split V = TLsh from 7i = Tig + Tih, and 
assume that the system-bath interaction is weak. That is, we will treat V = TCsh as 
a perturbation (we shall quantify later what we mean by "weak"; be patient!). We want 
to calculate U{t;to) perturbatively \n V = TisW, to that end, we will use the following 
Kubo-type identity [Ml p. 148]: 



■ i{w+y){t-to) ^ -iw(i-io) 







(3.11) 



This equatioE0 can be iterated by expressing e"^'/^)*^^^^)'^ in the integrand in just the 
same way, and so on. To second order in V one has: 

ft— to rt—to rt—to 



U{t;to) = e-^^(*"*°) 



• rt-to 1 /-t-to ft— to 

1 - - / ds V{s) - T2 / ds du V{u)V{s) 

^ JO "Jo Js 



(3.12) 



^ The decoherence approach does not exclude the existence of those superpositions. It only asserts 
that they vanish at too short a time to be observable [5B]. For a discussion of the state of the play in 
the search of quantum superpositions in macroscopic objects, see Leggett's insightful review |27j . On the 
other hand, decoherence does not seem to solve the "measurement problem". It is true that it provides a 
diagonal q, with its elements giving the outcome probabilities of each eigenstate. But it does not furnish the 
mechanism by which, is a given measurement, the system would "collapse" in the corresponding subspace. 

* The proof of (|3.11|l follows multiplying both sides by e"'"R'^''~'°' and differentiating with respect to t. 
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where V{s) is the Heisenberg evolution under 7i (interaction picture), that is, V{s) = 
Q+[i/h)'HsY^-{i/ti)'Hs^ Taking the Hermitian conjugate of (j3.12p we get a similar expression 
for U^{t;tQ). Both results can be inserted now in (jS.lOp . giving the formal time evolution 
of Qtotit) = Ugtotito)U^. 

As we seek for a differential equation, we can differentiate gtot with respect to t. 
Keeping only terms to order V^, one gets: 



^ = -^[n,gtot]-^[V,Qtot] + ^ 1^ *"ds (v[gtot,V{s-t)] + [V{s-t),gtot]V^ (3.13) 



where the operators without argument are assumed in the Schrodinger picture, whereas 
V{s — t) is the formal Heisenberg evolution under Ti.. The tilde in ^tot means unperturbed 
evolution under Ti, that is, Qtot = e"('/^)^(*~*o)^(to)e+('/^)^(*~*o). Keep in mind that 
^tot / Qtot- 

Now we have to trace over the bath to get the evolution equation for g, Eq. (I2.24|) . It 
is convenient to write the system-bath coupling as {F = F^, E = E^) 

nsb = F^E (3.14) 

with F acting only on the system's Hilbert space and E over the bath'sH Now we will 
assume that at the initial time to the density matrix Qtot is given by the following 
tensor product (decoupled initial conditions): 

Qtotito) = Qs Qh ■ (3.15) 

Besides, ^>b would be the equilibrium density matrix of the bath in the absence of the 
system, i.e., = e~^'^^/Z\y with Z\, = Trb(e~^^''). That is, the bath is considered at 
thermal equilibrium before "meeting" the system [see paragraph before (j2.1ip and (j2.12p ]. 
The factorized form ()3.15p is difficult of justify, but it is an assumption required to proceed 
further with the derivation. To be more reassured, think that in this work we are interested 
in the long time behavior, where the system would have "forgotten" any initial conditions. 

After all this preparations, let us trace over the bath's degrees of freedom, Trb(- • • ). 
Thanks to (j3.15p . we know that gtot{t) = gs{t) ® gh{t) because 7is and act on dif- 
ferent Hilbert spaces {[7{s,'Hb] = 0). Besides, g\^ oc dgh/dt = and gs{t) = 
^-ii/h)n4t-to) g^(t^^^Q+{i/K)n4t~to) _ Therefore, for the first order term of ([TTH]) we have 

^Trb i[V,gtotim = ^{EUF,Ut)] (3.16) 

where (• • • )b ^ Trb(- • • ^b)- In most cases of interest the bath is such that {E)\j = 0, and 
the above term drops. For the first term in the integrand of (I3.13P we have: 

Ti^iV[gtot,V{s-t)])= + Fg,{t)F{s-t)Tr^{EgkEis-t)) (3.17) 

- FF{s-t)gs{t)Ttk{EE{s-t)g^) 
= + Fg,{t)F{s-t)C{s-t)-FF{s-t)g,{t)C{t-s) 



^ The more general form Fj ® Ej is easily handled adding summations here and there. 
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where 

C(t) = {E{t)E)^ (3.18) 

is the correlator of the bath operator. This shows a clear analogy with the classical case 
[see Eq. (I2.14p ]. so we have used the same letter to denote it. 

Notice that to derive (|3.17p . the factorization (|3.15p has been essential; then Qtot{t) = 
gs{t) Qhit)- As the coupling Tish is factorized too, Eq. (|3.14p . we were able to take the 
parts acting on the system out of the trace action. Proceeding in the same way with the 
term [V{s — t), ^tot]^ in (|3.13p . we finally arrive at 

^ = -^[Ws,^>]+7^, (3.19) 

with the following relaxation term 

n = -^J^ ds^^^C{t-s)[F,F{s-t)Q\-C{s-t)[F,QF{s-t)]^ (3.20) 

Note that we have replaced by g inside TZ, as they differ in second order terms, while 
TZ is already of second order in V. 

Now we can be more specific on what we mean by weak interaction. The equation 
obtained, Eq. (j3.19p . was based on the expansion (j3.12p . There we assumed that V was 
a perturbation on 7i, so that it made sense to retain only the first few terms in the 
expansion. Well, but (j3.12p would only be valid at short times, Vt/h <^ 1, because letting 
1^ ~*o| ~^ oo the integral could become very large and break down the perturbation theory 
[34j . So, is it that the master equation only makes sense at short times? This would 
reduce considerably its applicability (let alone the description at short times being biased 
by the decoupled initial condition). Fortunately, the fact is that we implicitly assume the 
existence of correlation time Tc obeying 

C{tc) = 0, para t > . (3.21) 

Therefore, it is this Tc what limits the validity of the expansion as it bounds the integral 
term. Following Van Kampen [71 Chap. XVII] one introduces some 7 as a measure of 
the coupling strength (in the appropriate units [7] = [t""*^]), and asserts that the master 
equation (|3.19p is valid if 

7Tc<l, (3.22) 
quantifying in this way what we mean by weak-coupling. 

3.3 Heisenberg evolution 

In much the same way as we obtained an equation governing the evolution of the reduced 
density matrix, we can obtain the corresponding equation for any operator A within the 
Heisenberg picture. 

The Heisenberg evolution for A is given by A{t) = U~^{t,to)A{to)U{t,to)- We can 
use the same approximate expressions (j3.12p for U and and then differentiate with 
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respect to t. Keeping terms to second order in V, we would get an equation dA/dt = (•••) 
analogue to (|3.13p . Assuming then (i) the factorized initial conditions p.lSp . (ii) that 
the operator A only acts on the system's variables A = A(S'\ (lb is the identity in the 
Hilbert space of the bath), (iii) taking into account that A(t) = Tii,{A{t) q^) , and finally 
(iv) taking the partial trace Trb(- • • Qh), one arrives at 

d A i 

^ = -[n,,A]+n (3.23) 

with the following relaxation term 

n = -^ J^' ds\^C{s-t)F{s)[F,A]-C{t-s)[F,A]F{s)j . (3.24) 

Here the operators without time argument are understood to be evaluated at t. Note 
that TZ in the Heisenberg equation is structurally different (simpler if you want) from 
the TZ obtained for the density matrix g, Eq. ()3.20p . Remember that now we are in the 
Heisenberg picture, so that Qtot{t) = £)tot(0) = & ^b; we also exploited this to trace over 
the bath above. Finally, from the A solving the above equation one can get the desired 
average as 

{A) = Tt,{A0s) . (3.25) 



3.4 application to the bath-of-oscillators model 

Let us now particularize the master equation (j3.19p to the bath of oscillators model ()2.4p . 
Note first that the Hamiltonian ()2.4p can be rewritten as: 

Wtot= +F^^U^Qa + ^Y.^Pa+^lQl) (3.26) 

ns+^Fj:^\u^\yuji — ' ^ — . ' 

^ E Hi, 

where the counter-term ^-F" X^a l^aP/'^a been absorbed into 7is [see the discussion 
after (j2.4p and (j2.9p ] . In this form, the model has the structure Tish = F ®E for which we 
derived the master equation (I3.19p . (The case of Heisenberg evolution ()3.23|) is analogous.) 
In the quantum case it is convenient to rewrite (j3.26p as: 

Htot =ns + F®Y, (C"^^ + 4^-) + ^c^htba (3.27) 

a a. 

where h'^^ha are bosonic operators (creation and destruction) with the new constants 
Ca = \J /i/ (2cja) Uq,. We recognize the last term in (j3.27p as the standard quantum oscillator 
Hamiltonian tkjb'^ba- 

Let us now particularize the equations of the previous section. The correlator ()3.18p 
follows from the time evolution of the bosonic operators ba oc e'"^"* as 

a 
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where = l/(e^^°' — 1) is the Bose distribution. Introducing now the spectral den- 
sity (j2.13p as in the classical case: 



J{u) = vr - a;,) = vr J] ^5(0. - (3.28) 



we can write the correlator as 

C(r) = n — J{lo) [coth {^(3huj) cos(wt) - isin(a;T)] . (3.29) 
Jo ^ 

This renders the classical result (j2.14p in the limit 0, letting coth(/3/ia;/2) 2/{Phuj). 

Now we can write the master equation for the bosonic bath model ()2.4p to second order 
in perturbation theory 

^ = -^[w„^.]+7^ (3.30) 

with [cf. Eq. (fOn]) ] 

^ = -^ {i7(0)[F', 0] + I'ds C{t - s)[F,F{s - t)Q] - C{s - t)[F, gF{s - t)]| . (3.31) 

Here we have added to TZ the term coming from the counter-term in (I3.26p . with 7(0) = 
/i^^«Q/2a;^ which can be expressed in terms of the so-called damping kernel 7(r) = 
h Jq°° da;J(a;)/a;cos(u;r). That term cancels with i'y{0)[F'^ , g] arising after integration by 
parts in TZ. Therefore, as discussed in the classical case, all renormalizations on the system 
have been accounted for from the outset. 

Let us close with a few comments. Equation (I3.30p is sometimes criticized for not being 
of Lindblad type, and hence for not ensuring positivity in the time evolution of g. E In 
appendix D. 3 of PAPER IV we addressed this point from a physical rather than a rigorous 
point of view. Let us summarize our position. Equation (|3.3Up was derived under several 
approximations, the most important being 7 Tc <C 1; outside its validity range the equation 
would violate positivity and what not. Besides, the derivation assumed decoupled initial 
conditions at to = —00 [see Eq. (I3.15p ]. Therefore, taking arbitrary initial conditions at 
t = (which are the only cases of reported violations of positivity), it seems natural that 
the approach may produce nonsense. 

So then, can we only study the evolutions under those conditions? Yes and no. Imagine 
we want to study the evolution from an equilibrium state at t = of system plus bath. 
If the total system is in some sense "ergodic" it does not matter the initial condition we 
used at to = ~cxd (for instance, factorized, or other); at t = we will have our Boltzmann 



® A Lindblad type equation can be written as (h.c. = Hermitian conjugate) 

dg/dt = ^^[^s, £>] - ^ h„,m.{gLmL„ + L,nL„g - 2L„gLm^ + h.c. (3.32) 

n.m 

This is the more general time-local equation preserving the positivity of g. Its derivation is not based 
on any microscopic model, but in the theory of dynamical semigroups [2^1 Sec. 3.2], and is popular in 
mathematical-physics quarters. 
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Plot= Ps » Pb 



p(t = 0) 



► 



"preparing" the initial conditions 



dynamics studied 



t = -0O 



t = 



Figure 3.2: Scheme of the protocol discussed in the text to work consistently with the 
master equations. We assume that bath and system first met a long long time ago: long 
enough so as the decoupled initial conditions no longer matter. Then, from t = onwards 
we manipulate TLs with fields, etc., to study the response of the system. 

initial condition. From that time onwards, we can manipulate the system with fields, etc. 
(that is, modifying via 7is) and study the corresponding response (see Fig. 13. 2p . And all 
this because we are respecting a continuous process coming all the way from to = — oo. 
Therefore, if kept within its range of validity (weak coupling), the equation would preserve 
positivity, hermiticity, etc. to the approximation considered, d 

3.5 a few examples 

Within the oscillator bath model, there are some cases where we can write the master 
equation (I3.30p explicitly. 

3.5.1 when the eigenstates of Hs are known 

If we know the eigenvalues and eigenvectors |n) of the system Hamiltonian, TCs\n) = 
e„|n), we can write explicitly the equation for the components of the density matrix Qnm 
in that basis. To this end, we sandwich both sides of ()3.19p as {n\dtQ\m), and employ: 



where Anm = ^n — ^m- After some algebra (done with care in appendix D.l of PAPER IV), 
one obtains 



Now we can understand why if we use arbitrary initial conditions at t = with the constant coefficients 
of the asymptotic range {t > 0), we can get any kind of funny result. The g(to) allowed would only be 
those compatible with the original initial conditions at to — — oo plus the joint system-bath evolution. 

® Note that in the derivation we used a time independent Tis, while the manipulation protocol discussed 
is based on modifying Tis- There is no problem about this: as such time dependence results in a time 
dependence of some of the equation coefficients, but the basic structure remains the same [35) . 




(3.33) 




(3.34) 



n'm' 
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with the explicit relaxation term 

^ ^ T^nn'mm' Qn'm' ^ ^ [ ^mm' (X^i-^/|n''^™''^''^') 

""^ nm _|_ (r„|„/ + rj^i^,) F„/„Fmm' (3.35) 

— ^nn' {Yll^l\m'Pm'lFlm) ] Qn'm' ■ 

Here the relaxation rates are 

r„|^ := r(A„„) = i7(0) + / *° dTe-^^-^C(T) = VF(A„„) + i$(A„„) (3.36) 

Jo 

with the part 7(0) = Yla ''^a/^a coming from the counter-term [see after Eq. ()3.3ip ]. When 
the correlator C(t) is given by Eq. (|3.29|) . one can prove the important relations 

W{A) = e-^^T^(-A) , $(0) = (3.37) 

The first is the detailed balance condition. The second is a consequence of having accounted 
for the renormalization (and will be exploited later on). The equation obtained is of the 
Redfield type, like those introduced for the study of magnetic resonance [36j. 

3.5.2 RWA 

The "rotating wave approximation" is common in many branches of physics. Roughly 
speaking, the RWA consists in dropping the terms which "rotate" much faster than the 
characteristic frequencies of the problem addressed. As an illustration imagine we have 
terms like e'('^i+'^2)t ^nd e^^'^i-'^^)* (^^^h wi,u;2 > 0). Then one argues that eK'^i+^i)^ 
oscillates much faster than e'^'^i"'^^)* g^j^^^ g^j^ ]-,g replaced by its average (zero). Let us see 
how the RWA is implemented in our problems. 

Let us jump into the "rotating frame": "Qnmi^) — Qnm^^^^^'^"''^^ where the density- 
matrix equation (I3.34p reads: 

'^'^nmji) _ TP - ^(A„™-A„,„,)t ^^o oox 

— / ^ l^nn'mm' Qn'm'^ • \O.OOj 

n'm' 

The crudest approximation assumes that the terms e^l^^^^^^^^'m')^ with A^^ — ^n'm' 7^ 
revolve in time much faster than the characteristic relaxation scales of the system, so they 
are averaged to zero. That is, we discard all terms with A„m — ^n'm' 7^ 0. If we have a 
system spectrum without degenerations or regularities (see below), that is, where all A's 
are different, we have A„m — ^n'm' = 0; only if n = n' & m = m' or if n = m & n' = m' 
[55] . Therefore, in such a RWA the master equation splits into diagonal and off-diagonal 
elements as follows 

^Qnn 



n' 

<^Qnm i 



^ '^ T^nn'nn' Qn'n' (3.39) 
^nmQnm ~l" T^nnmmQn 



at n 



where we have returned to the unrotated basis. 
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It is not difficult to see that TZnn'nn' is real: TZnn'nn' = 2|F„„/pVF(An„'). Besides, 
T^nnmm = T^mmnn (^^ Q IS Hermitiaii). But in the RWA one only retains the real parts 
of TZnnmm, arguing that the imaginary parts reflect an extra renormalization. Notice that 
the equations for the off-diagonal terms are now uncoupled, and their solution is simply 

Qnm{t) = Qnm{0) e(-t'^-+^™)* (3.40) 

Therefore, we just need to solve the equation for the diagonal terms. Then analytical ad- 
vances are sometimes possible, and in the worst case, it will be computationally "cheaper" . 
Furthermore, equation ()3.39p is of Lindblad type [28] (by some chance), and the equilib- 
rium solution dtQ = results to be the canonical equilibrium oc e"^"^" (see next chapter), 
what makes it more "user friendly". 

As for our discussion of decoherence, given that the W{Anm) > (see PAPER IV), 
the TZnnmm give directly the inverses of the decoherence times. In this case the basis 
where the density matrix becomes diagonal, the pointer basis, is that of eigenstates of Tig 
(this is related with the equilibrium solution being oc e~^'^'' ) . Finally, the equation for the 
diagonal terms has the same form as the Pauli master equation, but there the off-diagonal 
dynamics is neglected. 

The RWA is many times difficult to justify (for a illustration of the limitations see 
[37|). For instance, we could have systems with Anm — A^/^^' very small. Or it could also 
happen that A„m — An'^m' = in more cases than those considered, for example, if there 
is a degeneracy, or if Tig has an equispaced spectrum Ann±i = cte. In the latter case, we 
would have A„m — An\m' = also when n — m = n' — ml for any shift p in the indexes. 
Then, we are not allowed to drop the Tlnn+p^mm+p terms in (|3.34p and we write 



d-Qnm i 

dt ~ ~h 



AnmQnm ~\~ ^ ^ T^nn+p,mm+p Qn+p,m+p (^•^-'-) 



We will refer to this equation as "improved RWA", and it will be invoked in the last 
chapter. 



3.5.3 semiclassical or high temperature equation 

This is another instance where the relaxation term can be written in a simple form. Let 
us start having a look at the formula (j3.29p for the correlator. Although the frequency 
integral extends to infinity, we have J{uj > i^c) = with ujc the bath cutofflf] Then, at 
high enough temperatures, Phujc ^ 1 one can approximate coth(Phuj/2) = 2/{Pfiuj). 

Now let us consider a low frequency part J{uj) ~ mjuj, which corresponds to the 
Markovian case discussed in the classical limit (see Sec. I2.3.ip . Under these conditions, 

^ Remember that in the classical case we wrote J{uj) ~ uj°' (see the footnote in Sec. 12.3.1) . But this is 
an idealized situation; in fact any physical J{u}) should go to zero at high enough frequencies |T| Sec. 3.1]. 
Therefore, one inserts a regularizator, exponential for example 

J{lo) oc u'^e-"^^" (3.42) 

and at frequencies higher than the cutoff the spectral density goes to zero. 
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the correlator ()2.14p is left as 



C(t) = 2m7 kBT5{T) + im% -^6{t) (3.43) 

dr 

Therefore the relaxation term reduces to 

^ = -^^^^1^ {[F,Fe] +h.c.) + g ([F, [ns,F]0] +h.c.) (3.44) 

This is the semiclassical or high-temperature equation. It has the advantage of not having 
an integral in TZ (time local) and hence is easier to handle. 

If we now particularize the above TZ to the case of a mechanical Hamiltonian Tig = 
/2m + V{x), with the particle linearly coupled to the bath coordinates, F = —x, we get 

^ = - —^^—[^^ e]] - i^I^' Ip^ eh] ' (3-45) 

with [A, B]^ = AB-\-BA the anticommutator. This is nothing but the "famous" Caldeira- 
Leggett equation, originally derived with the path integral formalism |9j. 

Let us close with a few remarks on the validity range as compared to the original 
equation (|3.3U|) . First, we used J{id) = rwyuj, with 7 having units of and being 
related with the coupling strength [see Eq. (12.130 ]. On the other hand, for such J{lo) the 
bath correlation time is Tc = Ph (see \28\ Sec. 3.6]) and hence the validity range from 
Eq. (I3.22P would be Ph'j <C 1. But on top of this restriction, Caldeira-Leggett put the 
PhuJc < 1 condition. Then, one fears that a priori the validity range would be quite 
restrictive. However, in the case of the harmonic oscillator. Tig = /2m + ^mJl^x'^, exact 
quantum master equations have been obtained [38l EO]. And for this case Unruh and 
Zurek [39j checked empirically that, except at short times, the restrictions on the high- 
temperature equation could effectively be relaxed to ujc/^ <^ Q^^T/h'y ^ Although limited 
to the harmonic oscillator, this is anyway good news on the applicability range of the 
Caldeira-Leggett master equation. 



3.6 summary 

We are aware that this chapter has been terse. Let us summarize the main points. 

First we have discussed the concept of quantum decoherence. The decoherence program 
is an example of how quantum theory can save face. In a way similar to how reversible 
dynamics could account for the irreversibility found in nature, by invoking the interaction 
with the environment, quantum mechanics manages to explain the "non-observation" of 
certain odd superpositions in the macroscopic world (for the good of cats and other crea- 
tures). The reason, we are told, is the extremely fast lost of coherence when we only look 
at our system, through the notion of partial tracing [Eq. (j2.24|) ]. 

Then we discussed the description of the dynamics of open systems by means of a 
master equation for the reduced density matrix q. Those master equations are valid in 
the range 7 Tc ^ 1 with 7 quantifying the system-bath coupling strength and Tc the bath 
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correlation time. The prime difficulty with such equations is in dealing with the time 
integral in the relaxation term TZ [Eq. (j3.20p ]. We saw how this can be handled without 
problems when the spectrum of Tig is known, getting a Redfield-type equation [Eq. ()3.35p ]. 

Finally, we discussed a couple of approximations widely used in the specialized lit- 
erature, like the rotating wave approximation [Eq. ()3.39p ] and the high-temperature ap- 
proximation a la Caldeira-Leggett [Eq. ()3.45p ]. Both equations have validity ranges, in 
principle, more restrictive than <C 1, although when they can be used result in a great 
simplification. 
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Chapter 4 

equilibrium properties 



Up to this point we were not concerned with the equihbrium properties in open systems. 
The equations discussed, classical or quantal, describe an irreversible dynamics in which 
the system in contact with the reservoir exchanges energy (chapters [2] and [3]) . Waiting 
long enough the system will equilibrate to the bath, reaching a stationary regime. Which 
is the state reached is the question that we will try to answer in this chapter. 

We will start revisiting the classical case. Then we introduce the thermodynamic 
perturbation theory formalism, and expand the partition function in powers of the system- 
bath coupling, getting explicitly the partition function Z to second order. After this, we 
will check that the master equations derived in the previous chapter are thermodynamically 
consistent. That is, that they have as stationary solution the equilibrium distribution 
obtained here to the same order in the bath coupling. We conclude discussing how to 
compute equilibrium quantities in open systems, and illustrate this with the example of 
the harmonic oscillator. 



4.1 the classical limit 

The best way to address the equilibrium properties starting from the dynamics is to resort 
to the Fokker-Planck equations. The necessary condition for Weq to be a stationary 
solution is i9tW4q = ^ciW^q = 0, where £ci is a compact writing of the Fokker-Planck 
differential operator. If we focus on the particle case [Eq. (j2.16p ]. we can easily check that 
Uq-p^^ = 0. Furthermore, it happens that any initial condition, in the hmit t ^ oo, 
tends to the same solution [14J. Therefore, the system relaxes to the equilibrium Gibbs- 
Boltzmann distribution VFeq oc e~^^=. The dynamics of classical open systems places us 
in a comfortable theoretical framework where the Fokker-Planck equations capture the 
approach to equilibrium. 

^ The first two terms in (|2.16p . arising from the Poisson bracket, give zero when appiied to e~^'^=. 
Letting the other two terms (diffusion and reiaxation) act on e~''^° one finds they cancel by using 
{^lm)dp\ v exp (-/?pV2m)] = -ap[exp(-/3pV2m)]. For spin probiems £ciWcq ~ foliows anaiogousiy 
from Eq. (fT^ . see 
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Figure 4.1: Pictorial representation of the change of variables Qa — UaxjuJ^ Qa where, 
for simphcity, we have used F = —x [see after Eq. (|4.3|) ]. Physically, what we do is to shift 
everything by a distance oc x. As the integrals extend from — oo to +00 their hmits do 
not change. Compare with the analogy of the Bohr- Van Leeuwen particle in a magnetic 
field (final comment in Sec. I4.2.T]) . 



Let us reexamine this from another point of view. Although in an effective/reduced 
way, the Fokker-Planck or Langevin equations describe the total dynamics of system 
plus bath (j2.4p (see the discussion on irreversibility in Sec. 12. 4p . But the total system 
is a macroscopic entity, and hence we can apply the Gibbsian hypothesis. That is, at 
sufficiently long times the averages can be computed with the Gibbs distribution for the 
whole system 



Q—PHtot 



tot 



^tot = J dpdx J YldPaYldQae-"'^'"' (4.1) 



-2^tot 

As we are discussing the reduced distribution, we integrate out the bath, getting 

f -i—r -I— r Q—fiHtot Q—fSHs 

We^= II dP, H dQ^— = (4.2) 



or in other words, Ztot = with = / H dPa H dQofi ^"^^^ and Zs = f dpdxe' 

This result has relied on the integrals over the bath being Gaussian (see Fig. 14. ip 

Ztot J 

so that the change of variables Qa + UaF/Lo^ — > Qa gives Weq = e~^^=/2s- This is the 
reason behing ()4.2p giving the reduced distribution equal to that of the system without 
couping. Now we understand that e~^'^''/Zs is stationary solution of ()2.16p . since the 
Fokker-Planck equation is nothing but the evolution of the total system with the bath 
variables integrated out. 

We remark that, although it results natural to accept that the system relaxes to the 
distribution e~^'^^, from the point of view of ()4.2p this is a "happy" coincidence. We 
mean that the result is "model dependent" and for other models (coupling nonlinear in 
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the bath coordinates, or other kinds of baths), the reduced distribution could be different 
from e~^^=. Otherwise, every system with several degrees of freedom would factor into 
the product of reduced distributions for each degree of freedom. 

In fact, this result is so particular that, as we will soon see, does not hold for the 
same Hamiltonian in the quantum case, and the system relaxes to a distribution different 
from the canonical one. What is the meaning behind this? In the books of statistical 
mechanics one assumes that the body is macroscopic and that the interaction with the 
bath is somehow mediated by its "surface" . As the surface L?' to volume ratio goes to 
zero in large systems oc the interaction with the reservoir Ji^h is put aside, leaving the 
Gibbs distribution [40j. But for the open systems we are treating here we cannot neglect 
TLsh- System plus bath, on the other hand, constitute a macroscopic object. And here it 
is where we can assume Gibbs. But then nothing ensures that the reduced distribution is 
going to be oc e~^^=. We insist again, the result ()4.2p is quite particular. 



4.2 the quantum case 

We now address the quantum problem. We will try to answer to the following questions 

• Is the reduced distribution the canonical one for the system without coupling? 

^eq := Trb k 6"^^= (4.4) 

To which we have already advanced the answer (previous paragraph), and, if it is 
not so, how is that distribution? 

• Is the reduced density matrix a stationary solution of the master equations discussed 
in chapter [3p That is 

A &q = . (4.5) 

with vCq a compact way of writing the master operator in ()3.30p . 



4.2.1 thermodynamic perturbation theory 

In principle we do not have a practical expression for e"^'^*"* (much as in the dynamics 
we did not know an exact differential equation for the reduced density matrix). Therefore, 
we will assume again that the system-bath coupling is weak. This will allow us to perform 
the equilibrium calculation to the same perturbative order as we did in chapter [3] for the 
master equation. Besides, with this calculation we will gain some insight in the way the 
first corrections to the uncoupled distribution emerge. 

In order to obtain an approximate expression for e~^^*°* we will start again from the 
Kubo identity (jS.lip . As in Sec. 13.21 we will identify 7i = 7is + Tih and TL^h = V and 
expand to the second order: 



-m 



1 



daV{-iha)+ / da / dO V {-iha)V {-ihO) 



(4.6) 
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This is nothing but the evolution (I3.12|) in imaginary time, with V{—iha) = e^'^^VeT''^'' . 
But note that here the expansion is in powers of f3V . 

We now trace with respect to the reservoir and apply the result to the bosonic bath 
model (as in section [3^ . getting 

Se, = '— |^e-/^^= + R (4.7) 




Here = Trs(e-^^=) and the "rest" R 

rP rcT 1-13 



R 



Zs 



[ da [ d9F{-iha)F{-ih{a -9))C{-ih9) + j{0) [ daF^i-iha) , 
Jo Jo Jo J 

(4.8) 

involves the correlator (13.180 in imaginary time C{—ia) = {E{—ia)E)\y. We already defined 
7(0) = ^Ea«a/2w^ after (fOHj) . which is a byproduct of the counter-term in ()2.4p . 

The result above is already properly normalized, Tr(£)) = 1, and Z^'^'^ is essentially the 



second order correction to Z^ot = Tr(e 
%^ = Zs+z(2) ; =-l/J^|F_|V^- [rdae'^^-C(-i^)-7(0) 



(4.9) 



The superindexes (2) tell that this result is valid to second order. The first order terms 
disappear as they did in the dynamics since in the bosonic bath (j3.27|) we have {E)\^ = 0. 

We have thus answered the first question, showing that the reduced distribution, 
contrary to the classical counterpart, does not agree with the canonical result at zero 
coupling. 

We close this section with a few more remarks on the result obtained. First notice that 

-A 



j d^7e'^^C(-a)-7(0) = cI>(-A) z(2) = -i/3^|F„„|V'3^"ci>(- 



nm) 



(4.10) 

That is, the equilibrium quantities, as computed from the partition function, will depend 
on the bath through the imaginary parts of the coefficients TZnn'mm' [Eq- (I3.36p ]. The real 
parts are associated to relaxation and decoherence times (section l6.2p . that is, to how fast 
we reach the stationary solution. But now we are in the statics, and we do not care about 
the speed of the approach, but about the final state instead; it then makes sense that the 
equilibrium results do not depend on the real parts of Tlnn'^mm'- 

It is instructive to see that in the limit fih — > 0, both Z^"^^ and the term R go to 
zero, recovering the classical result of damping-independent equilibrium properties in a 
Gaussian bath. Consistently, in the high temperature limit the quantum corrections to 
the canonical result disappear too (-Ztot = ZsZ\y). H 



^ In the high temperature hmit, this is understandable, since letting T — ^ oo we have -Ztot Ntot, with 
A'tot the number of states. But A'^tot = N^Nh so that Ztot = Z^Zh- Physically, at very high T the coupling 
potentials become negligible compared to the kinetic energies. 
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We conclude with an analogy to a famous problem: the absence of orbital magnetism 
in classical statistical mechanics. Let us start from the Hamiltonian of a charged particle 
in a magnetic field 

n = ^(p-eAf + V{q) 
Zm \ J 

and compute the partition function Z = j dqf dpe~^^. As the momentum integrals are 
Gaussian, on making the variable change Pi = Pi — eAi, the vector potential disappears 
from the final result, producing a zero magnetization oc dlnZ/dBz = 0. This is the 
content of the Bohr- Van Leeuwen theorem. 

In our case we have Gaussian integrals too, but it is Ua what disappears from Z (in- 
stead of A), so that the equilibrium quantities become independent of the system-bath 
couplings. Quantum mechanically (we are thinking in the obtainment of Z from the den- 
sity matrix), the Gaussian argument does not hold. Correspondingly, the equilibrium 
quantities do depend on TCsh, much as the charged particle developed a non-zero magne- 
tization (Landau). 

4.2.2 stationary solution of the master equation 

Let us address now the second question, that is, let us see if the damping-dependent equi- 
librium distribution ()4.7p is a stationary solution of the quantum master equation ()3.30p . 
We would not hide that the answer is in the affirmative. It should be so if we did things 
properly, since [e~^^'°' , Titot] = 0, while the master equations are just the Von Neumann 
evolution of the total system with some tracing over the bath degrees of freedom (remem- 
ber the classical discussion in Sect. I4.l"j) . What we aim is to prove here is that the equality 
still holds when e~^^*°* and Wtot are calculated order by order. Let us see thisi 

As all our calculations are perturbative (both in the dynamics and in equilibrium), we 

7 

will check CqQeq = order by order. Taking into account the decomposition (j4.7p of the 
density matrix g = + g^ we obtain the following consistency conditions 

(i) 

^ = ^[^^i?,Ws]=0 (4.11) 

(ii) 

^ = ^[^>g^Ws]+7^[4)] = o (4.12) 

Condition (i) is trivially fulfilled, being g^ a function of Tis- In order to prove (ii), note 
first that [Us, {Z^'^^/Z^) e"'^^^] = 0, so that we need to check only if 

[ns,R]+n[gl^] = (4.13) 

On writing R = Rnm\n-) {"ml we have [Hs,R] = Anmrnmln) {m\ with A„m the energy 
difference between the n and m eigenstates of TCg- Therefore, the equality (14.130 follows 



See [41] for a closely related calculation. 
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for the diagonal elements (n = m), given that 7^[iJ„„]|n)(n| = (we are omitting the 
marker "eq" in g to alleviate notation). Finally, is not difficult to prove the relation 

^i^nm = -^^(n|7^[^(°)]|m) (4.14) 

showing that the second condition is obeyed too. 

So we have just managed to prove the thermodynamic consistency of our master equa- 
tions. This is also relevant for the linear- algebra eigenvalue problem associated to the 
master equation. Solving it, one of the eigenvalues will be zero, corresponding to the 
stationary solution (|4.7|) . as we have just checked. As for the rest, we expect they having 
negative real parts, ensuring convergence at long times to the stationary solution. Positive 
real parts would imply a divergence with time which would not be physically acceptable. 

Well, at this point we have already answered our two questions. That is, the stationary 
solution of the master equation ()3.30p is the equilibrium result ()4.7p . and this does not 
agree in general with the canonical distribution at zero coupling. In other words, although 
by now it feels quite natural, the equilibrium properties of quantum open systems do 
depend on the damping strength. When we started this chapter this would have struck 
us, due to the intuition from the classical case and the Gibbsian prejudice. Indeed, in 
many works it is still stated that the system would approach the canonical e~^'^''. Here 
we have argued that in general this is not the case|l| 



4.3 computing thermodynamic quantities 

Once reached this point, we could use our formula (14. 9p for Z and compute thermodynamic 
quantities. To this end, it will prove convenient to redefine the free energy as follows (all 
other quantities will follow from this) 

^ = .Ftot-^b. (4.15) 

Here ^tot = —Tln{Ztot) with 2^tot = Tr(e~^''^'°') the partition function of the whole sys- 
tem, while .Fb = — Tln(2^b) with Z]^ would be the partition function of the bath uncoupled 
from the system Z]^ = Trb(e~^^''). That is, we are redefining the zero of free energy, re- 
moving the contribution of the bare bath. Besides, in this way our definition above for 
Z = Ztot/Z\^ allows obtaining this by the custom rules 

T = -kBTlnZ ■ Z = Ztot/^b ■ (4.16) 

The idea will be to use now for Z the perturbative result ()4.9p . 

From the definition of J- we can derive energies and entropies in the usual way 

U = F -TOtJ"; S = -dTT. (4.17) 

* In the original Lindblad calculation, the convergence to e~'"^° is used as one of the constrains to fix 
the form of the master equation from "first principles" . 
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Thus, our convention gives U = ^tot — and S = 5tot — since all quantities follow from 
T by linear operations. Then, in processes not involving the bath, we are just redefining 
the zero of all thermodynamical quantities. But even in the bath is involved, as we have 
merely extracted its bare contribution, it can be easily restored to the final results. 

Let us illustrate what these definitions would give in a calculation of the specific heat C. 
Imagine we want to find the C of our system, which as usual in this work, is an open system 
and interacts, for example, with the phonons of the lattice where it "lives" . To that end 
one applies a controlled heat pulse with a resistor, and measures the temperature change 
Ctot = AQ/AT. To isolate the contribution of the system from other sources (phonons), 
one can repeat the measurement in the absence of "system" (for instance, if it consists of 
magnetic ions, we measure a compound without spins, but having the same lattice). This 
would give the bath specific heat Cbj which experimentalists routinely substract from the 
original C. But the quantity C = Ctot — Ch is just what directly comes from our redefinition 
of the energy, when computed in the standard way C = —OtU- 

Instead of getting the thermodynamical energy as in (j4.17p . we could have assumed 
that the relevant system energy is {Tis) = Tr^gT-Cs). But in general this would not agree 
with our Z// 7^ (^s)- H The bare average of the system Hamiltonian 7is misses the important 
contribution of the system-bath coupling. And we do not mean, as above, the additive 
part of the bath alone, but the work done to "insert" the system in the bath. Indeed, the 
use of (TCs) to discuss thermodynamical principles instead of U can produce "violations" 
of the second law (see the lucid discussion of Ford et al. [l2l HSj Bil [35] ). 

On the other hand, for observables not depending on Tish — we have in mind the 
magnetization of a paramagnet where the coupling does not depend on the magnetic 
field — we have 



That is, in this case both the definition from the free energy, or directly from the reduced 
density matrix g, provide us with two equivalent ways of computing the magnetization 
(that we will use in chapter [8]). 

4.4 example: the harmonic oscillator 

Let us apply the thermodynamic formalism discussed to the example of a harmonic os- 
cillator Tig = p'^/2m + ^mfl'^x'^, bilinearly coupled to the bath F = —x. This exactly 
solvable model has been thoroughly studied in the literature [U [TSl US] . That allows us 
to compare our perturbative results with exact ones. As F = —x can be expressed as 
oc (a"*" + a) and the harmonic oscillator's spectrum is equispaced Ann±i = =F^^) we have 




(4.18) 



^(2) 



h 



(4.19) 



2mQ 



n 



n 



^ Such a way of computing the energy is possibly inherited from the "reductionistic" view of using the 
partial trace so extended in the field of open systems. 
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with en = hi}{n + 1/2). Doing the geometric sums we arrive at 



+ coth (Iphn) $+ (4.20) 



with $± = ^{hn) ± <^{-hn) and the $ function from ([336]) . 

Now we can obtain directly via (|4.16|) . and by derivation U, S, C, etc, as in (|4.17|) . To 
ihustrate, we obtain the following specific heat in the high-temperature range Ml/k-QT <C 1 

^ = 1 - 4^ (4.21) 

k-B 2-KkBT ^ ' 

while in the regime of low temperatures hQ/k^T ^ 1: 



kB~ KkBTj +3 nn^ ^^-^^^ 

Both results agree with the expansion to first order in the damping of the exact results 
by Ingold and Hanggi [37] . [fl 

Notice that in the low temperature regime, the dissipative corrections give C oc T. 
This is is contrast with the celebrated 7 = result where C drops exponentially at low 
temperature, due to the gap between the ground state and the first excited state. To 
make sense of this, we repeat again that now the total system is our system plus the bath, 
which removes the gap (the reservoir has infinitely many oscillators with a quasi-continuous 
distribution of Wq's). 

Another way of looking at this is resorting to the definition of density of states 

fc r+ioo+c 

d{e) = — d/3Z(/?)e^^ (4.24) 

27ri J-ioo+c 

That is, the inverse Laplace transform of the partition function. For the central harmonic 
oscillator, Hanke and Zwerger jl9] obtained d{e) exactly, finding that between the ground 
state and the first excited state there exists a finite density of states (that is, there is no 
gap, see Fig. I4.2p 



® We have used an Ohmic bath spectral distribution J(cij) ~ 77170; regularized with an exponential 
cutoff ujc- Then we know the explicit form of $: 

<1>(A) = 7 - 

in terms of the function [35] and the bath cutoff luc- 



2 In 



/ A 



V 27rfcBT 



1 + 1 + In — 

2-KkBT J \ A 



(4.23) 
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Figure 4.2: Density of states (j4.24p for the harmonic oscillator (figure extracted from 
[TS]). Only the ground state contributes a delta to (i(e), the other get "blurred" by the 
coupling with the bath, a feature often used when reasoning about open systems. 

4.5 final comments 

Performing the RWA in the master equation one gets Eq. ()3.39p . whose stationary solution 
coincides with the canonicallll That is, the equilibrium distribution at zero coupling. 
Possibly this result has contributed to the "belief" that the system always equilibrates 
towards e"^^'' in the quantum case as well. But in general this is not so, and here we have 
studied the actual form of the equilibrium distribution, and its relation with the stationary 
solutions of the quantum master equations. 

In the low damping range, the corrections computed here to the thermodynamical 
quantities, though sensible, are probably small and not easy to detect experimentally. In 
a way, this contrasts with the situation in dynamics, where weak coupling can produce 
sizable effects via relaxation and decoherence. Anyway, from the conceptual point of view 
the rigorous treatment of the equilibrium properties contributes to the understanding of 
those dissipative contributions, absent in the classical case. 



This can be seen noting that the non-diagonal terms go to zero. Inserting e = in the diagonal part, 
and using detailed balance p.37p one sees that e~'^'^° is indeed a stationary solution. 
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Chapter 5 

working in phase space 



We conclude the general discussion of the formalism of open systems with the reformulation 
of the master equations for q in phase space. The phase space approach, together with 
canonical quantization and path integral are three equivalent ways of quantizing [50]. In 
this formalism the theory "lives" in the space defined by the canonical conjugate variables 
of Hamiltonian dynamics (phase-space). Given that quantum and classical theory share 
now the same mathematical space, the formalism provides an ideal framework to study 
the transition from quantum to classical (or the other way around). In this frame we will 
be able to relate explicitly the classical and quantum theories of dissipation (chapters [2] 
and ED. 

We will motivate this chapter as Wigner did, in 1932, in his seminal article [51]. We 
will introduce historically the development of the theory for spinless particles. Afterwards 
we will provide a more modern view, following the postulates of Stratonovich to map 
Hilbert operators and functions in phase space. We then address the problem of systems 
with spin. In both cases (particles and spins) we transform the master equations into 
phase space. Then the corresponding Fokker-Planck equations are obtained in the limit 

"n^o". 

Keeping in mind the spirit of Caldeira and Leggett, who posed the theory of dissipation 
as a quantization problem [10], this chapter closes the circle, getting now the classical limit 
from the quantum case. The formalism discussed here will be very useful in Chap. [71 where 
we will solve quantum master equations in phase space by adapting techniques developed 
to solve classical Fokker-Planck equations. 

5.1 the problem 

In classical mechanics the dynamics is given by the Liouville equation: 

dtW = {n, W} 

where W and 7i are functions of the point v in some phase-space V. To illustrate, v = {x, p) 
in the one-dimensional particle case (V = M?) while v = [6, (j)) for classical "spins" (V = 
with §^ the sphere in three dimensions). 
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In quantum mechanics, the Liouville equation is substituted by Von Neumann equa- 
tion: 

where now TC and g are Hermitian operators in some Hilbert space. 

In spite of the formal similarities { , } [ , ] these two equations seem difficult to 
relate. We mean, it is not straightforward to take some appropriate limit in Von Neumann 
equation and arrive at Liouville's. They live in different mathematical spaces! In this sense 
it turns to be involved to study the classical-to-quantum transition, or the quantum-to- 
classical, in a continuous way. 

To bypass this problem we resort to the phase-space formalism, in which both classical 
and quantum mechanics "live" in the same space. Formulating the quantum theory as 
closely as possible to the classical case would allow to study their crossover, and also see 
how the corrections to the classical case are built; as well as studying the problems from 
another point of view. In particular, one can export notions and tools of the classical case 
to the quantum world. 

5.2 a constructive tale 

We are going to review the construction of the phase-space formulation of quantum me- 
chanics from a historical perspective. 

5.2.1 Wigner's function (1932) 

In that year Wigner introduces the following "distribution" [ST] : 



1 

W^{x,p) = ^J_ dye-'yP/'^rix - iy)V'(x + ly) , (5.1) 

where ip is the quantum-mechanical wave function. We have added the subindex to in 
order to distinguish it from the classical probability distributions (note that = W^). 
Next introducing position "eigenstates" \x), one defines the transform of any operator A 
as follows 

/oo 
dye™/^x-iy|yl|x + iy) (5.2) 
-oo 

Note that Wa_ = ^^+) that is, Wa is real for Hermitean A = A'^. With these definitions 
the expectation value of A in the state \^) can be written as an "average" over phase 
space: 

{iP\A\i;) = [dp [ dxWA{x,p)W4x,p) (5.3) 



Thus we see that plays the role of a distribution. Advancing results to fix ideas, 
Wa{x,p) is the classical dynamical variable A{x,p) corresponding to the operator A. 
Putting A = 1 (the identity), we have the normalization J dp J dxW^{x,p) = 1. 
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The transformation properties of are as follows [52]. Applying ■iIj{x) iIj{x + r) 
we have W^{x,p) — > W,^{x + r,p). When ip^x) ip{—x), transforms as W^{x,p) — > 
W^{—x, —p). Finally 'ip{x) i^*{x) gives W^{x,p) W^{x, —p). 

We thus see that is real, normalized, respects the symmetries, and the averages 
are computed like in classical mechanics (same for VF^,)!!] But the analogy is not fully 
complete, some differences should show up! is not positive, and besides < vr/^. 
That is, the value at any point {x,p) of the quantum distribution is bounded. This 
precludes, for instance, that one could write W^{x^p) = 8{x—xi)5{p—pi) for a distribution 
localized at (xi,pi). As Moyal pointed out in 1949, this reflects the uncertainty principle, 
forbidding well defined x and p "simultaneously" [53] . 



5.2.2 the 40s: Groenewold & Moyal 

In 1946 Groenewold expresses the transform of the product of two operators 

/oo 
^py/f^i^x-\y\BA\x + \y)iiy (5.6) 
-oo 

as follows 

op ox ox Op 

where the arrows indicate on which W the derivatives act upon. 

In modern terminology the transform of the product defines an inner product in V, 
called star product: 

Wba = Wb*Wa (5.8) 

Then the algebra of the product of operators (composition) is transformed in the ^-product 
in V, inheriting its properties. Namely, it will be associative, Wc(ba) = W{cb)Aj but in 
general non-commutative, Wab / Wba- 

Von Neumann equation gets transformed in phase-space as 

^ = -^[H,q] ^ dtW, = -^{WH*W,-W,*WH) (5.9) 

Introducing now Tis = p^/2m + V{x) and using Groenewold's result (15.71) we obtain Von 
Neumann dynamics in phase space 

dtW{x, p,t)=\-^d^ + v'dp + V w{x, p, t) . (5.10) 

L m ^ (2s -I- 1)! ^ J 

3 = 1 ^ ' 



^ One can proceed analogously with q (mixed states) introducing the analogue to (|5.1[1 : 

^^(■^'P) = ^J dye'^^^''g{x-y,x+ly) (5.4) 
So the calculation of averages (|5.3p now reads: 

Tr{AQ) = [dp f dxWA{x,p)Wg{x,p) . (5.5) 
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To abbreviate we removed the subindex q for the density matrix transform, and written 
y{2s+i) ^ Q^s+iy^ Besides, we have extracted the s = term from the sum {y'dp) 
gathering ah ^-dependent terms. In this way the formal hmit — > recovers the classical 
Poisson bracket: 

dtW{x,p,t) = {-p/md^ + V'dp)W{x,p,t) = {n,W] . (5.11) 

With the glasses of phase-space, the equation for the quantum dynamics is built from the 
classical one simply adding more terms to the Poisson bracket H 

Note that for the harmonic oscillator V = we have 1/(2*+^) = when s > 1, 

so Von Neumann's equation in phase space is identical to the Liouville equation. What 
is going on? Quantum and classical mechanics are the same for the harmonic oscillator? 
No, what happens is that in the classical and quantum cases the admissible solutions are 
different. Recall that W could take negative values and in addition is bounded [55J. 

In 1949 Moyal rounds up the work noting that Wigner's distribution (j5.2|) is just the 
inverse of Weyl's prescription for symmetric quantization ordering |4QJ, that is: 

xp ]^{xp + px) ^y^"" W(xv^'px)l2 = xp . (5.12) 

Besides, he computes the uncertainty relations in phase space |53j. Therefore, reached 
this point — calculation of observables, ^-product, and dynamics — the formalism in phase 
space was completed, constituting an alternative to the standard canonical quantization 
in Hilbert space. 

From the 50s onwards many applications were developed of this formalism. Many of 
them exploited its analogy with the classical formalism, adapting its schemes and tools to 
the quantum world j52l [561I571I58] . 

5.2.3 Bopp (1961) 

In 1961 Bopp introduces an operationally better way of computing the star product 

' ft d ' fij 

Wb*Wa = B{Q,V)Wa; ^-=^ + 75^' ^-=^"^5^ ^^-^^^ 

where B{Q^V) acts on the functions in V, with the same functional form a B{x,p) but 
with X p being replaced by the differential operators Q and V. As an aside, note that 
these operators obey in V the same commutation rules as their Hilbert analogues x and 
p, that is, [Q,V] = ih. | 

Let us apply Bopp's result to Von Neumann dynamics. To this end we use g = g'^ 
and TC = 7i^, whence (TCg)^ = gTL, together with Wa = for adjointing. Then, 

^ In this sense, the *-product can be seen as a "deformation" of the Poisson bracket. This motivated 
the inverse problem, that is, given a Poisson bracket how should we deform it to get quantum mechanics, 
deformated quantization; see [54] for an introduction to this program. 

^ This can be checked directly from the definitions (|5.13p of the Bopp's or by simply noting that 
W^[4,p1a = ihWA. 



52 



5.3. OTHER SYMMETRIES AND STRATONOVICH'S POSTULATES 



^He — ^qH^ ^^'^ obtain 

dtW = 2lui{rL{Q, P)) W . (5.14) 

The Hamiltonian 7i is now a function of the Bopp operators Q and V which act on the 
functions in . The equation above is a compact way of writing the overwhelming ()5.10p . 



5.2.4 Caldeira— Leggett equation in phase space (1983) 

Let us apply the formalism to the dynamics of open systems. This will prove useful in 
Chapter [71 We want to transform into phase space the master equation ()3.45p , which we 
briefly recall 

= --^[^s, q] p — [x, [x, Q\] - i—[x, [p, e]+\ , 

with [ , ]+ the anticommutator. To this end, the Bopp operators (j5.13|) are very convenient, 
since 

[x, [x, g]] = x^Q - 2xgx + gx^ > — > - c.c.j W = "^^^^ (5-15) 
(c.c. = complex conjugate) and 

[x,[p,g]+] = xpg+xgp—pgx — gpx \ — > (^QV —c.c?jW +[q,V* — c.c?jW = \h—pW (5.16) 

Introducing these two results in the Caldeira-Leggett equation, and accounting for the 
transformation (j5.10p of the unitary part —(\/h)\Hs^ g\^ one arrives at 

dtW{x,p,t) = \-^d:, + V'dp+^dp{p+mk^Tdj)+f^ ^f''^'^l[^ ^^^ . 

s=\ 

(5.17) 

We have written the equation is such a way that the first three terms correspond to the 
classical Klein-Kramers equation ()2.16p . And the last term contains the purely quantum 
terms coming from the closed evolution. It is important to remark that the terms (j5.15p 
and (|5.16p are the diffusion and dissipation terms in the classical case. That is, the 
Caldeira-Leggett equation is nothing but the Von Neumann evolution augmented with 
the classical irreversible terms. This identification has been possible thanks to the phase 
space formalism. Who could have guessed that [x, [x, ^)]] and [x, [p, were indeed the 
dissipative terms of the classical Klein-Kramers equation? 



5.3 other symmetries and Stratonovich's postulates 

We have seen how the phase-space formalism developed from the introduction of the 
Wigner distribution. That is, everything followed from the map between operators and 
dynamical variables in phase space [Eqs. ()5.2p and (15. 4p ]. But the map is restricted in the 
sense of "only" being applicable to "particle systems", we mean, to the algebra \x^p\ = i/i. 
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corresponding to the Poisson bracket {x,p} = 1. Then, how could we proceed if we want to 
tackle problems with other commutation rules? We have in mind a phase-space description 
for systems with angular-momentum type algebra [Si, Sj] = iheijkSk, corresponding to the 
classical Poisson bracket {mi,mj} = geij^Tn/^ of Eq. (j2.18p {g is the gyromagnetic ratio, 
and eijk the completely anti-symmetric unit tensor). 

The answer was given by Stratonovich [59j, completed by Varilly & Gracia-Bondia 
[_6(J], and later generalized by Briff and Mann for any finite Lie group In brief, the 

idea is to find a map between Hilbert operators and phase-space functions by imposing 
the "physically" reasonable properties obeyed by the Wigner function (section f5.2.1|) . 

5.3.1 the Stratonovich postulates (1956) 

Mathematically one starts with a physical system having a given "symmetry group" 
{[x,p] = ih, or [Si, Sj] = iheijkSk, etc) corresponding to a classical phase space V. □ Then 
one seeks for a map between the Hilbert space and V obeying the following Stratonovich 
postulates: 

(i) linearity 

A 1-^ Wa is linear and bijective. 

(ii) "reality" 

Wa+ = (Wa)* where yl+ is the adjoint of A. 

(iii) Normalization (d^ is the area element in V): 

/ WA{v)dfiiv) = Tt{A). 
Jv 

(iv) trace property 

/ WAiv)WBiv)dniv) = Tr{AB). 
Jv 

(v) group covariance 

WA4v) = WA{g-v) , ygeG. 

where G is our symmetry group and A^ := lA{g^^)AlA{g) within the group generators 
in Hilbert space. 

In [601 [61] it was proved that these properties determine uniquely the map. The map 
is bijective (i) and hence invertible. Condition (ii) ensures that Wa is real ii A = A'^ 
[see after (|5.2|) ]. Postulate (iii) is simply the normalization condition [see below (j5.3|) ]. 
Linearity plus (iv) ensure that the averages are computed as in Eq. ()5.3p . Finally (v) 
expresses respect to the basic symmetries of the system. In the case of particles, these 
postulates produce the Wigner function discussed above. 

* A somewhat more mathematical summary can be read in PAPER VI. 
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5.4 the case of spin systems 

Let us apply the formalism to spin problems. Here the dynamics is governed by a spin 
Hamiltonian Ti.{Sx, Sy, Sz) with 



[Si,Sj] — icijkSk 



(5.18) 



In the standard basis Sz\m) = m\m) and 5^|m) = S{S + l)|m), for m = —S, ■ ■ ■ , S, the 
spin quantum numberO In this case, following Ref. [60] we have the map 



WAie,^)=TriA{e,(l))A); A{e, 



Air 



2S I 



2S 



yE E TimYiU9,cl>), (5.19) 



1=0 m=-l 



where Tim = \/ (2Z + 1)/{2S + 1) x X^j S", i'l*?, j')(5, j| are irreducible tensors, 

the {S, j; S, j'\S, j') are Clebsch-Gordan coefficients and Y/m spherical harmonics [MIIIH]- 
Klimov and Espinoza [62] obtained a differential form for the spin ^-product which 
reads: 



(5.20) 



Here Ns = ^2S + 1 and is the angular momentum operator on the sphere 



d"^ d_ 1 d'^ 



C'Yim = l{l + l)Yi 



Im 



(5.21) 



Besides 



GiC^)Ylm = G{l)Yim , 



G{1) = V(2S + / + 1)1(25-0! , 



(5.22) 



fc=0 



do sin 6* ( 



(-ly 



j!(25 + j + l) 



(5.23) 



In formula (j5.20p the operators and G inside the brackets act only inside the bracket, 
while G"^ on the left acts on everything to its right. 

Equations (15.190 and (I5.20p are, respectively, the spin analogues of the Wigner map 
and Groenewold ^-product for particles. 



^ Recall that S = {fiB/g)/li with g the gyromagnetic ratio and /ib Bohr's magnetron (section I2.3.ip . 
In this section and in chapter |8] we will set hb ~ h = 1 (as in PAPER VI) . Then the classical limit will 
correspond to 5* — > oo. 
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5.4.1 Bopp operators for spins 



At first sight the ★-product ()5.20p looks of little use in explicit calculations. Then, in the 
spirit of Bopp, we would like to find a simpler representation for Wb * Wa- To this end 
we explicitly calculated in PAPER VI the transformation Ws^a, i = x,y,z, finding 



Ws,*WA=SiWA, 



x,y,z 



where 



Si = mir/i(£2) + i(i:n x L)if]2{C^) + \L, 



(5.24) 
(5.25) 



Here fji{C'^)Yim = i]i{l')Yim for i = 1,2, and the numerical coefficients are 



mil) 
mil) 



1 



2{2l + 1) L 
1 



(/ + 1)V(25 + 1)2 - (/ + 1)2 + 1^{2S + 1)2 + P 



V(25 + 1)2 - (/ + 1)2 - ^(25 + 1)2 + P 



2(2/ + 1 

The vector m{9,(j)) = (sin cos sin ^ sin cos 0) was introduced in ()2.18p while L is 

d 



(5.26) 
(5.27) 



1 m X 



d: 



m 



(5.28) 



obeying [Lj, Lj\ = ieij^Lj.. |f| Finally, invoking the associativity of the star product, for an 
arbitrary operator B{Sx, Sy, Sz) we have 



WBA = B{Sx,Sy,Sz)WA 



(5.30) 



where B{Sx,Sy,Sz) is the same function of its arguments as B{Sx, Sy, Sz), but now the 
Si are "classical" differential operators on the sphere. To abbreviate, instead of writing 
B{Sx,Sy,Sz) we simply write B{S). 

Thus the Si in the sphere play the same role as Q = x + \ih/2)dp and V = p — 1(^/2)5^ 
in M2_ Por example, in analogy with the Q and V commutation, we have here [iSi,^^] = 
iCijkSk- That is, the Si obey the same commutation rules as their Si partners (for the same 
reason as discussed in Sec. I5.2.3p . To conclude the analogy with the particle problem, the 
dynamics can also be written in the following compact form [cf. Eq. (j5.14p ] 



dtW = 21m( H{S))W 



^ The following formulas are of assistance 

[mi,Lj] = ; [(m x V)i,Lj] = ie,jfc(m x L)fe. 



(5.31) 



(5.29) 
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5.4.2 an example 

Let us apply the formalism to the spin Hamiltonian studied in this work 

Ti. = —DS^ — BzSz 
The Poisson bracket ()2.18p produces in this case the differential equation: 

dtW = - [IDS cos e + Bz]^ 



(5.32) 



(5.33) 



In the quantum case we employ ()5.30p and ()5.24p getting: 



d,W 



2D [cos efi^ + sin 9 ) + 



dW 



(5.34) 



The first thing one notices is that for D = (isotropic spin) the differential equation is 
the same in the classical and quantum cases: 



dfW = -B 



dW 



(5.35) 



much as in the harmonic oscillator problem (section I5.2.2p . But as discussed there, the 
admissible solutions are different in the classical and quantum problems. Second, when 
D ^ we can expand f]i and f]2 for large S [Eqs. (|5.26p and (|5.27p ] , obtaining 



(5.36) 



That is 7/1 = 25 and 7^2 — as S — > oo, recovering the classical equation (|5.33p for 
the anisotropic spin. Thus using the expansion (I5.36P one could get the corresponding 
semiclassical evolution with 1/S correctionsJil 



5.4.3 application to spin master equations 

As we did in the particle case (section I5.2.4p let us apply the phase-space formalism to 
open systems. We will consider the semiclassical regime, and eventually take the limit 
S" — > oo. Therefore, it will suffice to use the semiclassical master equation (j3.44p 

^ = -i[Hs, q] - Ar([F, Fe] - ^[F, [Hs, F]0] + h.c.) , (5.37) 

which has been written with the conventions used in this chapter (recall footnote [5]) . 

''' In PAPER VI the dynamics for a general quadratic Hamiltonian was obtained. Besides, we took the 
classical limit for an arbitrary Tis, getting the classical Poisson bracket as in the example discussed here. 
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In order to transform the master equation to phase space, we will resort to our spin 
Bopp operators (|5.24p . We need to compute [cf. (j5.15p and (|5.16p ]: 

[F,Fg]+ h.c. ^ (f{S) -ccVw (5.38) 



as well as 

[F,[ns,F]g]+li.c. ^ (f(cS)-c.c.)([W,(5),F(5)]-c.c.)w (5.39) 
Then the master equation ()5.37p gets transformed into 

dtW = |2Im(Ws(5)) +4At[W(F(5)) - ^Im(F(5))Im([?^s(5), F(5)])] jlV (5.40) 

As we did for the Caldeira-Leggett particle case, we will consider bilinear coupling, 
that is F = ^ ctiSi, linear in the system variables. In contrast with the mechanical case, 
here [Tig; -P'ldepends always of TCg, so that we cannot write the dissipative terms in a 
generic way|f| Considering the simplest case 7is = '^BiSi (isotropic spin), we obtain: 



1 d 



m 



= • 1 (m X Bcff) - m X A m X ( Beff - T— ) + M x Beff \W (5.41) 



d 



dm 



Here Befj = —dTCs/dm is the effective field of the classical equations, while we got an 
extra term M = m(f/i — S) +i(m x L)f/2. However, noting that fji — S = 0(1/ S) the term 
M disappears in the classical limit, recovering the Fokker-Planck equation (|2.22p . If] 

In this problem, even in the simplest case, the diffusion and relaxation terms contain 
quantum corrections. Recall that in the particle problem the diffusion and relaxation 
terms where identical to the classical ones [Eq. (I5.17p ]. 

There exists an analogy with classical spin problems. There a bilinear coupling induces 
multiplicative noise in the Langevin equation, in contrast with the particle case, where the 
noise enters additively for bilinear couplingl^ Formally, this is due to the angular Poisson 
bracket {mi,mj} = geijkmk which, even for F linear in mj, gives a non-constant {j4,F} 
in the Hamilton equation ()2.6p for any dynamical variable A |15j . Quantum mechanically, 
it is the angular-momentum commutation relations which render [F, ["Hg; -^]] non constant, 
which eventually produces the term M in the master equation. 



5.5 other distributions 

What we have told here is not the whole story. Canonical quantization comes associated 
with the ordering problem. In path integrals, similarly, it shows up in the choice of the 

Recall that in the generic case TLs = p^/2m + V{x) with F = — s, one has \Hb,F] = ihp/m, which 
does not depend on V{x). 

^ In PAPER VI we took the classical limit for generic Hb and F, also arriving at the Fokker-Planck 
equation (|2.22p . 

Compare the particle Langevin equation (|2.15p with the spin Langevin equation (I2.20p . In the former, 
the noise f{t) enters additively, whereas in the latter f{t) appears multiplying m and dF/dm. 
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evaluation point [63j. Analogously, in phase space this results in the ambiguity in the 
map between Hilbert operators and V functions — we can construct other distributions, 
each associated with a given ordering. In the 60s, Cahill and Glauber [6l] constructed 
distributions (called P and Q functions) in terms of the creation and destruction operators 
of the harmonic oscillator. Those distributions differ from the Wigner distribution (15. 1|) 
in that they correspond to normal and anti-normal ordering, respectively |52] . For normal 
ordering (P-function) [cf. Eq. (|5.12|) ]: 

^{xp + px) ^ P(xp+px)/2 = xp- 2ih (5.42) 

while for the Q-function P(xp+px)/2 = xp + 2ih. That is, the phase-space function xp 
corresponds to the operator xp in normal ordering and to px for anti-normal. 

These distributions can also be generated from the Stratonovich postulates of Sec 
15.3.11 Then the map gets parametrized by a superindex a and the transformed function 
is denoted by W^^; for o" = ±1 one has the normal and anti-normal cases, while o" = 
corresponds to the Weyl symmetric ordering discussed throughout. Postulate (iv), the 
trace property, is then modified as 

/ W^/\v)wif''\v)dfi{v) =Tt{AB). 
Jv 

Therefore the symmetric function (a = 0) is special, unique, in the sense that to compute 
the averages we only resort to functions with a = 0. Any other a requires transformed 
functions with a and —a. 

In the time before Varilly & Gracia-Bondia, the spin distributions followed the ap- 
proach of Cahill and Glauber. The construction was done, instead of using a and a"*", in 
terms of Bloch states or spin coherent states [65\ (those were the 70s). With this approach 
Takahashi and Shibata found the spin Bopp operators corresponding to the orderings 
a = ±1 |66j : an important result they called the product theorem. In PAPER VI, un- 
der the modern postulated-based approach, we generalized the formalism to any ordering 
(here we have described the a = case only). 

5.6 summary 

In this chapter we have surveyed the construction of the phase-space approach to quantum 
mechanics. The particle Bopp operators, and our generalization to spin problems, have 
allowed us to transform easily the master equations into phase space (sections 15.2.41 and 
I5.4.3p . With the equations so transformed we have recovered, in the classical limit, the 
corresponding Fokker-Planck equations, and we have seen how the first quantum correc- 
tions are generated. In this way we have connected the classical and quantum theories of 
dissipation of chapters E] and O On the other hand, as we will see in Chap. [71 transforming 
the Caldeira-Leggett equation (15.170 to phase space with be instrumental for its numerical 
solution. 
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Chapter 6 

methods 



This chapter is a transition between the general treatment hitherto discussed and the 
chapters to come, devoted to specific apphcations. Here we explain what we intend to 
calculate and how to carry it out practically. 

We present first the linear response theory for the obtainment of dynamical suscepti- 
bilities. Then we use the example of the two-level system to illustrate how those sus- 
ceptibilities can give us relaxation and decoherence times. Finally, we introduce the 
continued-fraction method to solve master equations and the practical computation of 
response functions. 

6.1 linear response theory (general formalism) 

We start recalling the protocol we discussed in chapter [3] for the study of the dynamics 
using master equations (section 13.41 fig. 13. 2p . During — oo < f < we let the system 
thermalize to the bath. At t = we add some field, or change some parameter in the 
system Hamiltonian Tig, and study the dynamics. Thus master equations are a natural 
frame to address problems where the system is perturbed, and study its adjustment to 
the new conditions. 

In order to not modifying the intrinsic nature of the system studied, the probe can be 
made small enough jG^. But then the dynamics can be studied within perturbation theory. 
Linear response theory (LRT), along with dealing with such adjustments, provides the link 
with a typical experimental situation — the measurement of the dynamical susceptibilities. 

So hands at work. For — oo < t < we can write the dynamics formally as 



where Cq is a compact way of writing the master equation with TCg = TCq for t < 0. |1| When 
t = the previous dynamics would have brought the system to the stationary solution 

^ Throughout this chapter we use the notation corresponding to the dynamics of g. The treatment, 
however, is general and hence vahd for the examples of cleissical Fokker-Planck equations of chapter [2] or 
their quantum counterparts in phase space (chapter [SJ . 



dtQ = CqQ , 
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(recall chapter H]) 

CqQ(, = . (6.2) 

Then at t = we modify TLs as Us = 'Ho + 5 ■ TCi, with 5 accounting for the size of the 
perturbation. To first order in 6 we can write: 

£(t >0) = £o + <5- A , (6.3) 
and the formal solution of the master equation is: 



g{t) = Qo + 6-0i{t); gi{t) = f ds e^'-'^^^ Ci{s)qo (6.4) 

J — oo 

We now ask about the evolution of the average of some observable (^4) = TiiyAg) and 
introduce: 

AA(t) := {A){t) - {A)o (6.5) 

where (^)o = Tr(A^)o) is the average in the absence of perturbation. Then /S.A{t) measures 
how far we are from the reference unperturbed evolution. 

The generic form of the change is TLs — > Ho + 6 ■ Ti.if{t) with f{t) some function of 
time. Correspondingly 6 ■ Ci{t) ^ 5 • Ci f{t), whence 

/oo 
dsR{t-s)f{s) (6.6) 
'OO 



The response function R follows from (j6.4p in the form 

Therefore our task will be to calculate R{t), which is independent of f{t) within linear 
response theory. We will carry this out in two typical situations: the relaxation experiment, 
and the AC response. 



6.1.1 relaxation experiment 

Let us consider a step function perturbation /(t) = Q{—t), so that f{t) = 1 for t < while 
f{t) = when t > 0. Note that we have reversed the "set up", saying that in the interval 
— oo < t < the system was "perturbed" with a constant 6, and at t = the perturbation 
is removed (Fig. 16. 1|) . H 

In this case {A)o = {A){oo), since the "unperturbed" evolution takes place for t > 0, 
and we know that the system will eventually reach the stationary solution. Then (j6.5p is 
written as 

Ar A{t) := {A) (t) - {A) {^) (6.8) 

That is, ArA{t) describes how the equilibrium is approached as t — > cxd. Hence the name 
"relaxation experiment" and the subindex "r" used. 

^ This is just a matter of terminology, as we could have said that Hb is unperturbed in the interval 
— oo < t < 0, while at t = we included the perturbation. Either way the result is the same. However, 
the language used in this section seems more natural, as we will see. 
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"i ^ ! ► -^l 1 KaA 

t=-0O t = '"'^ t = -'>3 t = 

Figure 6.1: Left: scheme of the relaxation experiment. We let the system equilibrate in 
the presence of a weak force 6 and at t = we switch it off. Right: AC experiment with 
a time periodic probe 6 ■ sin(u;i) applied from t = on. 



6.1.2 AC response (frequency domain) 



According to the experimentalists, it is often easier to measure the response to an oscil- 
lating probe than to conduct a relaxation experiment. In is convenient to introduce then 
the Fournier transform of the response function R: 



dte-"^'i?(i) 



At i = we switch on an oscillating probe f{t) = e" 
follows as 

AAc^(t) = {Am - (A)(0) = 6 ■ x(a;)e-* 



(6.9) 



see Fig. 16. ip and the response 



(6.10) 



Besides linear response theory teaches us that x(^) of the AC experiment is related with 
the response AAr{t) in a relaxation experiment [68, App. C]: 



X 



cq 



1 — icj 



AAr{0) 

Here x^'^ = x('^ = 0) is the response to an infinitely slow field, so that 



(6.11) 



X 



oq 



AAr{0)/6 



(6.12) 



6.2 Bloch equations 

To fix ideas and give content to the formalism, let us consider the solvable problem of the 
two-state system. In the basis of eigenstates |1) and |2) we can always write 

ns = -lB,a, (6.13) 

where az is the Pauli z matrix and the energy levels are ei = —Bz/2 and €2 = Bz/2. The 
notation reflects that we cannot help thinking of a ^ spin with a Zeeman Hamiltonian. In 
this way we bear in mind a physical realization of the two-level model. 
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The above Tig is coupled to a bosonic bath as in the model ()3.27p . The dynamics wiU 
be given in general by the master equation (|3.34p : 

~ J^^nmQnm ~^ ^ ^ T^nn' mm' Qn' m' (6-14) 
n'm' 

As we set the problem ourselves, we make our life easier assuming that the rotating wave 
approximation can be made (section l3. 5. 2p . Then the above system of ODEs simplifies; in 
particular the off-diagonal elements decouple and evolve trivially. For the two-state prob- 
lem the remainder equations for the diagonal [Eq. (13.390 ] are called the Bloch equations 

m. 

dgii/dt = -P2\iQii + Pi\2Q22 (6.15) 
d£»22/dt = +P2\lQll - Pl\2Q22 (6.16) 

Here Pi\2 = W^i2|-^i2p) with the matrix element of the coupling F{a) and the rate W, 
providing detailed balance Pi\2 = P2\i [Eq. (|3.37p ]. 

As mentioned the non-diagonal elements simply obey {qi2 = Q21) 

dQi2/dt = +^BzQi2 -TdQi2 (6.17) 

d£'2i/dt = -^B^Q2i -TdQ2i (6.18) 
with Yd = — i?ii22 = —^2211 > 0. The solution ()3.40p then reads 

^>i2(t) = ^>i2(0)e(^^^/^-r'*)S (6.19) 

and its conjugate for Q2i{t). Therefore one introduces the time constant T2 = l/T^, 
measuring the speed at which the off-diagonal elements go to zero — the decoherence 
time. 

We are left with solving the Bloch equations for the diagonal elements. Using that 
1 = Tr(^>) = gii + ^>22 and introducing the population difference = gn — g22, we have 

M^{t) = M,{0) ■ e-^'-* + M,{oo) , (6.20) 

with the relaxation constant = {Pi\2 + -f2|i)/2- (The notation again tells that we have 
in mind the spin problem, with the magnetization Mz oc {(Jz)-) Following the literature 
one defines the relaxation time Ti = l/F^, which measures how fast the diagonal terms 
equilibrate to their stationary values. And these, as we pointed out in 14.51 agree with the 
Boltzmann distribution M^(oo) oc e"^^^ — e"^'^^. H 

This solvable example, plus LRT, will illustrate how one can obtain Ti and T2 from 
the susceptibility curves. Given that the susceptibility is an experimentally accessible 
quantity, this provides a way of measuring relaxation and decoherence times. 



^ Recall that under the RWA the "pointer basis" are the eigenstates of TLs (chapter [H]). 
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6.2.1 LRT and Ti 

Let us reconsider the relaxation experiment of section 16.1.11 in a two-level system. For 
t < the system was perturbed by an extra constant dBz'- 

rLs = -\{B, + 5B,)a, . (6.21) 

At t = we remove the probe, SB^ = 0, and follow the evolution. 

Immediately before t = the system would have reached the stationary state in the 
presence of B^ + SBz- Under the RWA, we repeat again, the stationary solution is the 
canonical one: £»ii(0) = el^iB,+5B.) ^ ^^^^q) = Q-PiB.+5B,) ^^^^q^ ^ ^^^^q) = 0. The 

same at t = oo, but now with 6Bz = 0. As 5Bz is suitably small we can expand to first 
order getting Mz{Q) = Mz{oo) + 6Bz ■ dsMz- All the quantities here are equilibrium ones; 
by definition ObM^ = x'^'^ is the equilibrium susceptibility, so we can write 

AM,{t) = SBz-xl''e-^^' (6.22) 

We can plug the above relaxation function into Eq. (j6.11|) . getting the dynamical 
susceptibility Xz{^) as@ 



Xz 1 + i^^i + iu> 



(6.24) 



We have written the result for Xzi^-^) in two equivalent ways found in the literature. These 
functional forms are called after Debye's seminal work on dielectric relaxation. Let us split 
into the real and imaginary parts 

which have been plotted in Fig. 16.21 In the zero frequency limit the real part tends to 
the equilibrium susceptibility. Besides, the maximum in the imaginary part is located 
at w = l/Ti. In an ac experiment, where one applies a sinusoidal field and records the 
response, one can directly obtain x^^ a-nd Ti. 

Another usual form of representation is the Cole-Cole plot, with the imaginary part 
plotted vs. the real part (Fig l6.2l right). For Debye relaxation one has Im = -^/(l — Re)Re, 
in evident notations; that is, a semicircle with radius 1/2 located at (0,1/2). Note the 
marked points: the zeroes of the imaginary part are at w = and oo, corresponding to 
1 (maximum) and (minimum) of the real part. The maximum of the imaginary part, 
Im = 1/2, is reached at a; = 1/Ti, taking half the value of the maximum real part. This 
plot is used as indication of non-Debye behavior if deviations from a semicircle are found. 



Here we use 

roc A 

dte-(^+'"^* = — ^ . (6.23) 
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Figure 6.2: Left: real part (solid line) and imaginary part (dotted) of the Debye 



curve (j6.24p vs. frequency. The arrows mark special points in the spectrum, namely, 
the peak location of the imaginary part, which gives Ti, and the equilibrium susceptibility 
approached in the zero frequency limit. Right: Cole-Cole representation (Im vs. Re). 



6.2.2 LRT and To 



Let us now apply the LRT formalism to a problem with a transverse probing field (see 
also PAPER IV) . At negative times the Hamiltonian of our two-level system is then 

1 



-B^a^ + 5B 



(6.26) 



where cTx is the corresponding Pauli matrix (we could use (jy instead). At t = we switch 
6Bx off. As before, for the stationary solutions we will use the canonical form. Then to 
first order in 5Bx we have Qii{Q) = ^?ii(oo) = e^^^, Q22{^) = Q22{oo) = e~^^% for the 
diagonal elements, while the coherences (off-diagonals) read Qi2{0) = £'2i(0) = dB^Mz/ Bz 
and £'12(00) = £121(00) = 0. 

Now we monitorize ax (the operator coupled to the perturbation) through its average 
Mx = Q12 — Q2i- This relaxes from its equilibrium value oc Xx^ = M^/Bz as follows 



AM,(t) = 5BxX. 



eq 



(6.27) 



g(iB-rrf)t _^ ^{-iB-rd)t 

From this relaxation curve we can get the AC susceptibility Xx{'^) [Eqs. (j6.1ip and (|6.23p ]: 

X? i(u; + B) + T,^ iico -B) + r,- ^^"^^^ 

This is different from a simple Debye form, and we are going to see that includes absorption. 
Let us extract the real and imaginary parts 



Xx 



+ + Bz)Bz ^ + (u; - Bz)Bz 



+ i 



+ BY + Tl 
{u; + BY + rj 



(a; 



BY + rl 



(6.29) 



+ 



BY + Tl 
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Figure 6.3: Real part (solid) and imaginary part (dotted) of the transverse susceptibility 
Xxi^) [Eq. (I6.29P ]. The parameters are = 5 and Ta = 1- The arrows mark the 
equilibrium susceptibility and the width of the absorption peaks. 



These are plotted in Fig. 16.31 where we see again that the equilibrium result is recovered 
as w — > 0. As for the imaginary part, it consists of two Lor entzian- type curves, centred 
around the level difference ±Bz and having width F^. In this case the perturbation did 
not commute with the base Hamiltonian ()6.26p . inducing transitions between its levels. 
Spectroscopists call T2 (= l/T^) the life-time of the levels; for us it is the decoherence time 
(the time to erase the off-diagonals). Letting F^ ^ (closed system) the absorption lines 
go over two Dirac deltas at ibi?^. Their broadening in an open system gives T2 directly. H 



6.3 beyond the 2-state model 

The above example was a simple case where we could do everything analytically. In 
general, we will have to deal with quantum master equations like ()3.30p . But we saw that 
if Tis has a discrete and finite spectrum, we can write a density-matrix equation for the 
Qnm as in (j6.14p . This is a linear system with N x N equations. Diagonalizing the system's 
matrix (numerically most likely), we obtain N x N eigenvalues and eigenvectors. Then, 
with the appropriate initial conditions we can write 

NxN 

Qnmit) = Qnmi^) + Yl cS^'^^^e"^'* , (6.30) 
1=2 

where we have split the contribution of the zero eigenvalue (giving the stationary solution) . 

With these Qnm{t) we can construct the relaxation curve lS.A{t) for any observable A 
of interest [Eq. (j6.8p ]. Plugging it into Eq. (|6.1ip . and doing the integrals with (j6.23p . we 

^ The "magnetic resonance" experiment can be done alternatively by fixing uj and varying (the level 
splitting), looking for the resonance. 
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(»/Ai Re[)Cz]/Zz 

Figure 6.4: Left: real and imaginary parts of xi'^) given by the sum of two Debye terms, 
with weights oc 0.6 (left) and 0.4 (right). The arrows mark the relaxation times and the 
equilibrium limit. Right: the Cole-Cole plot (cf. Fig. 16. 2p is approximately given by two 
semicircles, as the two relaxation times are separated by several orders of magnitude. The 
zeros of the imaginary part correspond to the real parts 1, 0.4 and (see left panel). The 
maxima Im = 0.2 and 0.3 correspond to Re = 0.2 (= 0.4/2) and Re = 0.7 (= 0.4 + 0.6/2). 
The arrows mark these special points. 



obtain the corresponding dynamic susceptibility 



^^('-) = E«^A4h^' (6-31) 



The Qi oc m ''i"'"*^ (""1^1"^) involve the matrix elements of A (PAPER HI). Real Aj 
(the Tr in the two-level system) contribute Debye-type relaxation terms to the suscepti- 
bility. On the other hand, complex eigenvalues (like iF^^ it in the 2-level system) would 
contribute with absorption profiles. 

Let us consider now an observable such that the of the complex eigenvalues are all 
zero (the observable does not see those modes). Then the total susceptibility is made up 
Debye's; an example with two of them is plotted in Fig. 16.41 Each peak on the imaginary 
part is located at the corresponding Aj. 

In a multiexponential case like this, one tries to characterize the response with a few 
effective parameters. The integral relaxation time Tint is defined as the area below the 
relaxation curve: 

rint = ^ dt^^; x(c^) = X'''(l - i^Tint + • • •) (6.32) 

This quantity governs the low frequency behavior of the susceptibility. On the other hand, 
the high-frequency properties are better characterized by the initial (t = 0) slope of the 
relaxation curve. This defines the effective time Tgf 

_i d fAA{t)\ 



[aa{o)) 



■ x(^) = X^" (6.33) 

t=0 ^'^ef 
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Both time constants coincide in a problem with a single exponential decay. We will see 
examples with several relaxation/absorption peaks in Chapters [7] and [51 

6.4 continued-fraction methods 

It is not always possible, or practical, to diagonalize our master equations (large number 
of levels N, continuous spectrum, etc. . . ). In some cases we can resort to the continued- 
fraction method [H]. This method is a relative of the solution of kinetic equations by 
expansion into complete sets (a la Grad) [69]. It has been fruitfully exploited to solve 
classical Fokker-Planck equations [HI \7U\ and quantum master equations [TU [721 [ZSl [Zl] • 
The ideas is to expand the unknown solution q of our dynamical equation dtQ = Cg in 
some basis {uj} (all this applies to the classical or quantum W) 

Q = Y^ C,u^ . (6.34) 

i 

Plugging this into dtQ = Cg one gets a system of coupled differential equations for the 
expansion coefficients, that is: 

I 

Ci = Qi^i^jCj (6.35) 
j=-i 

Here / is the range of the coupling. If / = 0, then Ci = Qi^iCi, and Ci{t) = Ci{0)e'^^-*^ 

solves the problem — we had diagonalized the master equation! 

When / 7^ 0, on the other hand, one proceeds taking the Laplace transform g[s] = 
dte~^*g{t) (it is after all a system of ODEs with constant coefficients). Using the 

transform of the derivative g[s] = sg[s] — g{0), the system of differential equations is 

finally transformed into a set of algebraic equations for the Ci, with recurrences between 

the indexes 

- C,(0) = {Qi^i - s)Ci + Qi,i+jCj . (6.36) 

Compared with (I6.35|) . the coefficient of the central term was simply augmented to {Qi^i — 
s), while — Cj(0) appears as an inhomogeneous term. 

6.4.1 solving recurrence relations: case 1=1 

Let us assume that, with a smart basis choice, we have ended up with a recurrence 
like (|6.36|) with 1 = 1 (first non-trivial case). Then we introduce some simplified notation: 
fi = Ci(0) for the source terms ("forcing"), and the constants Q~ = Qj.j-ii Qt = Qi,i+i 
and Qi = Qi^i — s. Then the recurrence (|6.36|) would take the "canonical" form: 

Qrc,_i + Q^Ci + Q+a+i = -fi . (6.37) 

At this point we could diagonalize the matrix associated to the system ()6.37p . Then, 
instead of diagonalizing C directly, we would have made an analytical effort leading to 
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a tri-diagonal matrix; the numerical work would be more efficient. But we can take an 
analytical step further. If we introduce the following ansatz 



Ci = Sid^i + ai (6.38) 

in equation ()6.37p , we easily get the following relations for the "ladder" coefficients Si and 
the shifts aj (associated to the inhomogeneous term) 

Qf „. _ + Qfo'i±i (6 39) 



Qi + Si±i Qi + Si±i 

Now we can explain the name "continued fraction" given to this method; iterating the 
denominators with the same formula one could write: 

S^ = (6.40) 

Q. + Qt ^'-^ 

This would provide a "closed form" for the coefficients Ci[s]. We could say that this is a 
semi-analytical technique, where we end up evaluating (numerically) continued fractions. 
Eventually, inverse Laplace transformation would give the time evolution of the coefficients 
Ci, and hence of the full distribution g. H 

6.4.2 Brinkman hierarchy: matrix recurrence relations 

The best way to learn how to convert a kinetic equation into a suitable recurrence relation 
is through an example. Let us consider the Klein-Kramers equation (j2.16p . which in 
appropriate dimless units reads 

dtW{x,p)= -pdr, + V' dp + -fdp{p + dp) W{x,p) . 

We follow Risken and begin expanding in a basis [cf. Eq. (|6.34p ] 

W{x,p) = w{x,p)'^Cn{x,t)tpn{p) (6.41) 

with the pre-factor w extracted for later convenience. Inserting this expansion in the 
Fokker-Planck equation we get the following "equations of motion" for the coefficients c„: 

Cn = Qn,n+mCn+m ] Qn,m = j '^P tpnC-tpm ] £=W~^Cw (6.42) 

Here C is the differential operator on the right-hand side of the Klein-Kramers equation, 
while C is some similarity transformed version of it: Cf = w~^C{w f). 



® This last step may be seen as a handicap of the method, as the inverse Laplace transformation can 
be quite unstable [48], requiring many s-points to attain accurate results. As we will see later, in the cases 
treated here (statics and stationary ac responses) we do not need to use Laplace transformation. 
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Following Risken we put w = e (p^/2+^(^))/2^ which results in 

w-^ dp {p + dp)w iJm = (d^ + - ^) V''" • (6-43) 

That is, the dissipative part acts formally as the Schrodinger Hamiltonian of a harmonic 
oscillator (on the p dependence). This cries for introducing "ladder" differential operators 
a = dp+ p/2 and a'^ = dp — p/2 obeying [a, a+] = 1 and giving 

w'^ dp {p + dp)w ipm = a^a ipm ■ (6.44) 

On the other hand, the closed part (Poisson bracket) gets transformed as 

w-^ {-p + V' dp) w^^ = i-D+a + D_a+) i>m (6.45) 

with the shifted derivatives D± = dx^ V /2. Finally, if we choose for the V'n the Hermite 
functions, ipn{p) oc Hn{p/V2), the term a^a becomes diagonal, giving Qnn, while a and 
a+ produce Qnn±i [Eq. ()6.42p ]. Putting all this together one gets the so-called Brickmann 
hierarchy [TSl [76] 



Cn = - {VriD_Cn-i + jncn + Vn + ID+Cn+i) , (6.46) 

having the form of a 3-term recurrence. 

However, as the D± still contain operators on x, we expand on a x-basis {tia(x)}, 
obtaining a recurrence relation, now without operators, but involving two indexes: c" = 
Qn'n+m^n+m ■ Order to couvert this into a one-index recurrence (the problem we know 
how to solve), we introduce the following vectors and matrices: 



Cn = Q-Cn-1 + Q C„ + Q+C„+i . (6.47) 



Here L is some truncation index for the x-basis, and the matrices Q are 

Qa,a+I3 = Qn]n'^'^ = -nj 5a,a+l3 

with x-matrix elements Aa/3 = J dx UaAu/^. 

Finally, we could Laplace transform as in ()6.36p . getting an algebraic recurrence 
like (|6.37|) . with the Q's and C"s given by the vectors C and matrices Q. The resulting re- 
currence could be solved like the scalar case of subsection 16.4. II The only difference is that 
the fraction bars in (|6.4U|) have to be understood as matrix inversions {A/B — > B^^A). 
In the statics no Laplace transform is required, and Risken solved the problem studying 
parameter ranges out of the reach of other techniques [14 
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6.4.3 the case / > 1 

What happens if the recurrence involves 5 or more terms? Then we build appropriate 
matrices and vectors (by mounting the ones above) and the recurrence can always be cast 
into the "canonical" 3-term form [13] • The method will be useful as long as I is not too 
large. We will solve problems where / = 1 or 2. 

6.5 stationary solutions, LRT & continued fractions 

Let us see how to apply the continued-fraction method to obtain stationary solutions and 
dynamical susceptibilities. The approach will be used later in chapters [7] and El 

6.5.1 obtaining stationary solutions 

When we want to obtain the stationary solution £f? = we do not need to resort to 
Laplace's transform. Indeed, plugging the expansion q = CiUi in Cq = [cf. Eq. (|6.2p ]. 
we get 



for the coefficients of the stationary solution. This is an algebraic recurrence relation 
like (I6.37p . with fi = 0, solvable simply as discussed in Sect. 16.41 

6.5.2 application to LRT 

Let us have another look at the AC experiment of section 16.1.21 Recall that from t = 
on we apply a sinusoidal probe f{t) = e"^*; then the evolution operator becomes C = 
Cq + 5 ■ e"^*>Ci, while the solution can be written as ^> = + <^ " e"^*f)i. Compare with 
Eqs. (16. 3p and (16. 4p : here we have explicitly extracted the t dependence e''^* from q and 
C. Therefore, both qi and Ci are time independent (of course, qq and Cq too). 

Using now q = ia;£»ie"^* and equating terms with and without oscillating factors we get 



Inserting the expansion qq = ^^CiUi in the first equation we get the recurrence ()6.49p . 
which we can solve by continued fractions, getting the Ci. Next, we insert the solution 
for go in the second equation and expand qi as in (j6.34p . finding a recurrence like (j6.49p . 
but now with fi ^ 0. The inhomogeneous part comes from CiQq, thus justifying the 
carrying along of the forcing term from the generic recurrence (|6.37p . In summary, with 
two recurrence relations we can get g to first order in 6. With this g we obtain the AC 
response of the averages of interest [Eq. (j6.10p ]. 



I 




(6.49) 



i=-i 



= Cogo 
iujgi = Cogi+Cigo 



(6.50) 
(6.51) 
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6.6 final comments 

In this chapter we have considered the system's response to first order in the probing 
field/force. To go beyond LRT, one would need the response to higher orders [cf. Eq. (16. 4p ] 

Q = 00 + 5 ■ Qi + 6'^ ■ Q2 -\ (6.52) 

We will see in chapter [8] how to proceed along this line. 

In the second part of the chapter we have briefly introduced the continued-fraction 
method to solve master equations. This is a semi-analytic technique, and for the cases we 
are interested in (stationary and first few order responses) we will see that it can be very 
efficient. Besides, it can be used when the spectrum is continuous (Chap. [7]). In addition, 
one gets the full g, not only the observables. 

It has some drawbacks too. It is very specific of the problem addressed, requiring a 
specific pre-analysis in each case (finding a good basis {uj} to expand g, obtainment of the 
required matrix elements, etc.). Besides, it can become numerically unstable in certain 
parameter ranges (typically very low damping and temperatures |14j). 
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Chapter 7 

applications I: 
brownian particle in 
a periodic potential 



In this chapter we apply all the preceding open-system formalism to the problem of a 
quantum particle in a periodic potential V{x) = V{x + L). That is, we study the Bloch 
problem taking into account the coupling to a dissipative bath. Besides the substrate po- 
tential, the system can be acted upon by a force —x • f{t), allowing the study of transport 
properties (the current). The external force breaks periodicity and the band picture, so 
that the problem becomes a delicate one (continuous unbounded spectrum, Bloch oscilla- 
tions, Wannier-Stark ladders [20] . . . ). On top of this the particle interacts with the bath, 
further complicating matters. 

The chapter begins with the solving of the Caldeira-Leggett equation (I3.45|) in phase 
space [Eq. (j5.17p ] by means of continued fractions. In the classical limit this equation, 
which goes over the Klein-Kramers equation ()2.16p . was solved in this way by Risken (sec- 
tion [632]) • Here we discuss the extension of that work to the quantum regime (PAPER I). 
In the paradigmatic example of a cosine potential V{x) = cos(x) we will see how the classi- 
cal and quantum distributions W{x,p) can differ. After this we address transport problems 
by adding a driving —x ■ f{t). We compute the current and the mobility in the cosine 
potential as well as in potentials lacking space-inversion symmetry: ratchet potentials of 
the type V{x) = sin(x) -|- r/2sin(2x). The phase-space approach will allow us to connect 
smoothly with the classical results. 

7.1 solving the Caldeira-Leggett master 
equation by continued fractions 

We start writing the Caldeira-Leggett equation in phase space as [cf. Eq. (|5.17p ]: 

oo 

dtW{x,p)= \-pd, + v'dp + jTdp{p + dp)+Y,^^''^y^^'^'\x)dl^'+''^]wix,p) . (7.1) 
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The first two terms are simply the classical Poisson bracket. The third one includes 
the dissipation and diffusion originating from the coupling to the bath. Being a high 
temperature semiclassical equation (Sect. I3.5.3p . those terms are equal to the classical 
ones. The sum in the last term incorporates the quantum mechanical "corrections" to the 
closed Hamiltonian dynamics (Moyal terms, section fS. 2. 2p . 

Equation ()7.ip was made dimless by means of a characteristic length xq (e.g., the 
period of the potential) and a characteristic energy Eq (the potential amplitude) . Besides 
the parameters have been thermally rescaled introducing ujt = {k-QT /rnxo)^^'^ . Thus the 
momentum p in (|7.1|) was made dimless and thermalized (the potential V —> V/k-QT too), 
while the other parameters read 

^^-JL. rS') - ^-^y ( V ■ (7 2) 

^^"^t' ~ {2s + 1)1 [u;tk) ' 27T-h- ^^'^^ 

Here ojq and 5*0 are the characteristic frequency and action [ojq = (E'o/'^^^o)^^^ 
5*0 = Eq/luq]. The temperature has been absorbed in the definition of jt and the Moyal 
coefficients k^^\ In this way the form (j7.2p shows that the quantum terms become less 
important the higher K and/or the temperature are. Indeed the parameter K tells how 
quantum the system is {K (x 1/h, Chap. [21 Sect. 12. 6p . Thus the limit K ^ oo recovers 
the Klein-Kramers equation from the Caldeira-Leggett master equation in phase space. 
From now on we will simplify further the notation by omitting the subindex T en 7. 

7.1.1 the classical limit (K 00) 

As we stated in chapter \5\ one advantage of the phase-space formalism is that by varying 
one parameter the master equations go over their classical Fokker-Planck counterparts. 
In our problem, letting i^T —> cxd in ()7.ip we recover the Klein-Kramers equation ()2.16p . in 
dimensionless form 

dtW{x,p)= [-pd, + v'dp + -fdp{p + dp)'\w{x,p) (7.4) 

In section 16.4.21 we transformed this equation into a recurrence form, the so called Brick- 
mann hierarchy (I6.46|) . Recall that what we did was to expand the distribution in Hermite 
polynomials for p while the x-basis was left unspecified, W(t,x,p) oc c"{t)ua{x)^^Jn{p)■ 
T'his led to a 2-index recurrence = Q^n+m'^n+m which was converted into a recurrence 
in one index by defining the vectors (|6.47p and the matrix coefficients (|6.48p 

C„ = Q-C„„i + QC„ + Q+C„+i . (7.5) 

In this form, the recurrence could be solved by continued fractions (Sect. 16. 4p . The x- 
basis choice is usually dictated by the space symmetries of the problem including boundary 
conditions. 



^ Note the relation between K and the thermal de Broglie wavelength. The definition AdB = Aa; = h/Ap 
with Ap — \/Mk^T (thermal uncertainty from equipartition) gives 



AdB=-^o^/-g,;| (7.3) 
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application to periodic potentials 

If the substrate potential is periodic, one Fourier expands as V'{x) = Ylq^q^^'^^ ■ Then it 
is natural to choose plane waves for the x basis 



Uaix) 



27r 



a 



0,±1,±2, ... 



(7.6) 



qa 



(7.7) 



To get matrix elements = J dx UaAup, we use orthogonality J^^ dxe^'^^e 
obtaining 

{dx)a,a+l3 = ; (^')a,o+/3 = ^ • 

From these we can get the matrix elements of Q, which have the structure 
[ia6a,a+i3 i V^]. Hencc, the number of diagonals with non- vanishing elements equals the 
number of harmonics of V (only one for the cosine potential) . 

Hitherto we always had in mind solving the recurrence (|7.5p in the index n. But 



OC 



we could also have chosen to transform the 2-index recurrence 
introducing the following vectors 



o,a+/3 a+f3 



n+m^n+m 



by 



/3=-/ 



(7.8) 



In this case / equals the number of harmonics in V, since Qa,a+i3 = for /? larger than 
that number. Both recurrences (in p or in x) can be employed at convenience as Risken 
already did in the classical case [H]. 



example: the cosine potential 



Let us consider the stationary solutions W{t — > oo;x,p) for a particle in a cosine 
potential, tilted by a constant force (see fig. 17. 

V{x) = -Vo cos(x) -x-F (7.9) 

To this end we set C„ = (or = 0) and solve the vector recurrences. With the c's so 
obtained we reconstruct W{x,p) oc X] Ua(x)V'n(p) [cf. Eq. ()6.4ip ]. which is plotted in 
figure [7^2] (still classical results). 

Let us follow Risken in the interpretation of these graphs. In the absence of ap- 
plied forces {F = 0) the stationary solution is the equilibrium Boltzmann distribution 
oc e~P In the opposite case, when the force is very large, the substrate potential 

matters only a little, and we can think of a damped particle dragged by a constant force. 
In this limit (p) = F/j. Therefore, from small to large forces the distribution shifts in p, 
from the Boltzmann case with (p) = towards the stationary solution oc e~2(P~^/'>')^. In 
other words, a transition from "locked" to "running" solutions as a function of the force 
(with bistability in between). 
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Figure 7.1: Potential profile V{x) = — cos{x) — Fx. In the Hamiltonian case (7 = 0) the 
particle can be "locked", oscillating in one well, or "running" with {p) ^ 0. These two 
cases are depicted by solid trajectories. When coupling the particle to a bath it experiences 
fluctuations and dissipation. These can lead to the eventual trapping of a running solution 
into one well, or to the thermal launching of a locked particle over the potential barrier 
(dashed line) possibly into a running solution. 




Figure 7.2: Risken's graphs: classical distributions for a particle in the tilted periodic 
potential (I7.9p . The distribution is periodic in x so we have just plotted two periods for 
illustration. The damping is 7 = 0.05 and the forces F = 0, 0.075, 0.15, and 0.2 (left to 
right, top to bottom). Note the {p) = F/'y ~ 4 for the largest force, and the bistability 
between running and locked solutions for intermediate forces. 
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7.1.2 the quantum regime {K < oo) 

When K is finite, we have to take into account the rest of the terms in the Moyal series 
in equation ()7.ip . Here is when we start to appreciate the advantages of the phase-space 
formaUsm: one part of the equation is already solved — Risken did it for us! Now we are 
left "only" with handling the last sum, i.e., k^'^ V^'^'+'^\x) d^p'^'^^W{x,p). To this 

end, we make use of the result 



w 



-^dl'+^ [w i^m) = (a+ + a)2^+ Vm, (7.10) 



where recall that a = dp + p/2 and a'^ = dp — p/2. Then, using the Hermite p-basis 
each term in the sum will couple n with n it (2s + 1). That is [77], the Moyal term will 
in principle give coupling to all indexes in the quantum generalization of the Brickmann 
hierarchy (16.460 

m 

The explicit form of the matrices Q can be found in PAPER I, sections 4 & 5. It is 
important to remark that for polynomial potentials = 0, after some s. Then 

the recurrence would be finite, providing the finite I we required in the previous chapter 
to solve by continued fractions. On the other hand, for non-polynomial potentials there 
exist non-vanishing derivatives at all orders, like in a cosine potential. We are back in 
the generic case of the quantum Brickmann hierarchy ()7.1ip with infinite coupling range, 
making continued fractions of no use. But we still have the x-recurrence, let us see if we 
are more lucky with it . . . 

the case of periodic potentials 

In this case we proceed expanding the x dependence in plane waves (j7.6p . The "new part" 
(the Moyal sum in powers of h) brings in the derivatives ^^'^^ = V^^'^^e^'^^. To get their 
contribution to the matrix elements A^p = J dx UaAup, we use the following result in the 
plane- wave basis: 

[y^'\x)U^^^ = V^'K (7.12) 

Therefore we have Qn'n+m — whenever /3 is larger than the number of harmonics in 
the potential. Then, using the vector recurrence (j7.8p in the a index, we have a finite 
coupling range, just as in the classical case. The matrices Qa,a+f3 will contain more stuff, 
/i-dependent terms, but what is important is that the index I is finite and the same as in 
the classical problem: equal to the number of harmonics. In this case we could not chose 
between x or p recurrence (of no use), but luckily one of them works. 

example: cosine potential 

As we did in the classical case, let us compute the stationary solutions of the master 
equation in a cosine potential. This potential has one harmonic (I = 1) yielding a 3-term 
recurrence relation. Setting Cq, = and solving, we get the distributions of Fig. 17. 3i Here 
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Figure 7.3: Wigner distributions in the deep quantum regime {K = 1, left, & K = 2, 
right). For the temperature T = 1 and damping 7 = 10~^ the distributions are close 
to the canonical equilibrium (see text). The white islands correspond to zones where W 
becomes negative (slightly: W ~ — 10~^ in the minima of the islands closer to p = 0). 



we have set a "very" quantum regime, with K = 1,2 (recall K oc Sq/K). In order not to 
leave the validity range hj/ksT <C 1 of the Caldeira-Leggett master equation (Sect. [33T3|) . 
we have used dimless 7 = 10~^ and T = 1 (ensuring 'y/TK <C 1). 

The first striking feature (already announced in chapterfU)) is the distribution becoming 
negative in some regions (the white "islands" in figure 17131) . Besides, compared with the 
classical graphs 17. 2[ the distribution is quite delocalized in x. Specially for K = 1; for 
K = 2 we start to recognize the first symptoms of localization. 

To interpret with the hands the distributions of figure 17.31 we recall first that in the 
range 7/T = 10~^ <^ 1 the stationary solutions should be well approximated by canonical 
distributions (chapter H]). Recalling also how we build the Wigner function from the 
density matrix [Eq. (j5.4p ]. the equilibrium Weq would be: 



1 



2TTh 



/OO r 
dye^^P/Vq(x - + iy) ; Qe^= dq e-^'^\q){q\ (7.13) 
-00 J 



This will be approached by our f — > 00 stationary solutions. The Bloch waves \q), for small 
K, can be constructed perturbatively ([78], section "potential treated as a perturbation") 



|g) = e'''"(l + ?7e'" + r/*e-^") ; 



rj 



K 
2^ 



(7.14) 



Thus, when the characteristic action is "small" compared with h (small K/2t:) 
particle behaves almost as a free particle (see the energy bands in figure [73]). 
Inserting eventually the above \q) in the equilibrium ()7.13p we obtain: 



the Bloch 



V 

2\n\- 



2|^H +e 



(p+ir 



cos(2j;) 



V 
2M- 



+ 2\r] sin(x) 



2|.7|2 



+ e 



(7.15) 
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This remarkable form captures qualitatively the features of the numerically obtained 
distributions. The first term (a Gaussian in momentum) "is" the free particle e^"?^ — > 
exp(— p^/2|?7p), while the second comes from shifted plane waves e^^''^^-*^. They are re- 
sponsible for the delocalization along the x axis. The last two terms in (j7.15p follow from 
the cross products between e'^'^"''^^^ & gi('?-i)a; g^^^ qIq^ ^ gi(?±i)a^^ respectively. Specifically, 
the third term incorporates a modulation due to the potential profile (as it is of order |r/|^ 
is less noticeable for K = 1). The fourth term, finally, produces negative regions. 

The above terms, and hence the negative islands, can be seen as an expression of the 
wave character of quantum mechanics. By this we mean that W generates the averages 
and hence it is made of squared wave functions ['tp*{x — u/2)ip{x + u/2), section [5.2.1) . 
Seen in this way, W is like an intensity (its Fourier transform) and the last terms in ()7.15p 
depict a typical interference pattern between sinusoidal waves. 



7.2 computation of observables: transport problems 

We did not mention it in this chapter, but we recall that once the distribution W{x,p) is 
reconstructed from the expansion coefficients c's, we can compute any observable following 
essentially the classical prescription [Eq. (15. 3p ] 



{A) = j dp j dxA{x,p)W{x,p) . 



The "classical observable" A(x,p) is actually the Wa transform of the operator A in 
Hilbert space, Eq. ()5.2p . 

In some cases, it is not required to construct W explicitly, as some observables of 
interest can be expressed directly in terms of expansion coefficients. For example, in 
transport problems the current (p) can be written as 



(P) 



(7.16) 

n \ ct / 



The coefficients la and -K^2n+i depend on the prefactor w in Eq. (|6.4ip . and on properties 
of the integrals of Hermite polynomials and plane waves [given explicitly in Eqs. (29) and 
(30) of PAPER I]. 



7.3 example I: quantum transport in a cosine potential 

Let us focus now on transport problems, and the effects of the dissipative bath. We 
consider first that the substrate potential is the simplest possible, the cosine, which is 
however tilted by a constant force —x ■ F (see figure [7!T]) . For this problem we solve the 
Caldeira-Leggett equation and compute the current in stationary conditions {p){t — > oo). 

Figure 17.41 shows the results for several values of the damping. The phenomenology is 
that studied by Chen and Lebowitz, which did a difficult path-integral calculation treating 
the potential perturbatively [79| 180] . The curves can be understood as follows. When the 
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Figure 7.4: Left: average velocity {p) as a function of the applied force, for various 
7 in the weak damping regime. Right: energy bands of the (untilted) cosine potential 
V{x) = — Vocos(x) for K = 1. 

force is small the particle moves on the band bottom, which for K = 1 has a near free 
particle form (figure [7^ right). Then {p) = F/7 as in a damped and dragged classical free 
particle. Increasing the force, we approach the band limit. There Bragg scattering reduces 
the group velocity oc dE/dq^ and hence {p) decreases (see any solid state book, e.g. [SI])- 
But increasing further the force, we go into the next band, progressively recovering the free 
particle behavior {p) = F/7. Figure [73] also shows how increasing the damping the curves 
approach this free particle current. This is understood in terms of the folk argument of 
bath as blurring the energy levels (recall figure 14. 2p , which in our case amounts to bridge 
the band gap, responsible for the decrease. 

7.4 example II: directional motion in a ratchet potential 

Up to this point we have treated the cosine potential. In this section we consider the 
minimal extension of adding a second harmonic 

V{x) = -Vb [sinx + (r/2)sin(23;)] (7.17) 

with r the ratchet parameter (r = brings us back to the one-harmonic case). This 
potential is asymmetric, as seen in Figs. 12.31 and 17.51 Looking at the potential profile we 
realize that the slope to the right is less steep (easy side) than to the left (hard side). 

An important property of ratchet potentials is that they can accommodate directional 
motion. By this we mean that the system can have a net current even if the average 
of the external forces is zero (on average, the force does not do work). One also speaks 
of rectifying zero- mean perturbations (generating a net response). How is this possible? 
First, if the system is out of equilibrium, the second law of thermodynamics does not 
hold. Besides, we need the potential being asymmetric, as otherwise there would not 



82 



7.4. EXAMPLE II: DIRECTIONAL MOTION IN A RATCHET POTENTIAL 



2 



1.5 - 



1 

0.5 — \ - /' \ r r 

- 

-0.5 : \ / t / r~ 

-1 - v 

-0.5 -0.25 0.25 q 0.5 

Figure 7.5: Ratchet potential profile ()7.17p for r = 0.44, and tlie associated band structure 
for K = 5, having only two bands below the barrier top. 























-F 






2 Ji/co 




t 









Figure 7.6: Scheme of the square-wave force profile applied in the problem of rectification 
in a ratchet potential. 

exist a privileged direction [82^ [55] . Then a system modelized by an asymmetric potential 
like ()7.17p . and coupled to a dissipative bath, can show such a phenomenology if a driving 
force X • f{t) keeps it out of equilibrium. 

For the numerical solution we are going to assume, for simplicity, that f{t) has a square 
wave profile, which changes between positive and negative values zizF with a period 27r/a; 
(see Fig. 17. 6p . Besides, we use the adiabatic trick: the period is long enough (compared 
with the rest of time scales of the problem), so that most of the time inside each step +F 
or —F the current is the stationary one at the corresponding force. Then the net current 
(averaged over time) would be: 

1 {P)r = l {P)+F + 1 {p)-F , (7.18) 

with {p)±F the corresponding stationary solutions!! 

Working in the stationary regime, we set Ca = in the recurrence ()6.47p and we are 
back in the case discussed in section 16.5.11 As for the continued-fraction solution, note 
that having two harmonics in the potential, the associated recurrences now have / = 2 as 
coupling range (section I6.4.3P . 



^ Note that for the cosine problem we have {p)+f ~ — (p)-f and hence no rectification {p)r = 0. 
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Figure 7.7: Rectified current in the classical limit for the ratchet potential ()2.27p driven 
by the square wave of figure 17. 6[ Results shown for various values of the damping at 
T = 1. The thicker line is the analytical result [T3] in the overdamped regime 7 — > oo. 
The vertical lines mark the two critical forces ()7.19p . 

7.4.1 directed motion in the classical limit 

Let us try to understand first the classical problem in order to address later the effects of 
adding quantum effects. The classical rectified current in shown in figure 17^ The main 
feature is that the net current is positive, meaning that the particle moves on average to 
the right, the easy way. A way of understanding this is introducing the critical forces for 
barrier disappearance, to the left and to the right: 

|F+|/yo = l-r; |F-|/yo = l + r. (7.19) 

Thus, the barrier felt by the particle when tilting the potential with +F is lower than that 
when we tilt with —F, so it is easier to move to the right. 

To understand the other features, we will think first on the T = limit; that is, only 
dissipation, no fluctuations. The associated deterministic system has two more charac- 
teristic forces Fi and F2, obeying Fi < F2 < F^. When F < Fi the dynamically stable 
solution (the attractor) is the locked solution, while for F > F2 the running solution is 
the stable one. Well, if the system is locked close to a minimum (F < Fi), or if it is in a 
wild running mode, in both cases the potential asymmetry (and other details) become less 
relevant. Therefore, we understand that the rectification will be more noticeable in some 
intermediate force range, some "window", as seen in figure rm Besides, the characteristic 
forces Fi and F2 are damping dependent, and happen to decrease with 7 |I3J. Thus the 
window of maximum rectification shifts down along the force axis when we decrease 7 |84j . 
Finally, to restore T in the picture, we can think of its effect as allowing jumping between 
attractors, smoothing the sharp ranges, and in particular allowing rectification below F^ 
(PAPER I and PAPER II). 
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Figure 7.8: Rectified current in the quantum problem, for several values of oc Sq/H. 
We fixed the damping to 7 = 0.05 and plotted the classical curve {K — > 00) as reference. 

7.4.2 quantum corrections 

Our goal here is to study how the rectified currents of the previous section are modified 
when the quantum properties of the system are taken into account. To this end we focus 
on one of the curves (we fix one 7) and progressively decrease the quantumness parameter 
K oc Sq/K. This makes quantum effects more and more important. 

In figure EH] we see that for large amplitudes of the force wave \F\ the quantum cor- 
rections increase the rectified current + {p)-F, whereas the rectification is decreased 
at low forces with respect to the classical case. 

To understand the reduction/increase, in PAPER I and PAPER II we resorted to the 
semiclassical calculation of the transmission through a saw-tooth potential. This is an 
asymmetric potential, much as the ratchet potential, with easy and hard sides (figure [7l9|) . 
but where analytical calculations are possible. We studied the transmission coefficient 
when deforming the potential with -\-F and —F (appendix F, PAPER I). We found 
that for energies over the barrier the transmission is higher to the easy side, because the 
wave-reflection is less intense than to the hard side (schematically indicated by the broken 
arrows in figure [72]) • On the other hand, when the particle has an energy lower than the 
barrier, the transmission is higher to the hard side, as tunneling is more probable towards 
that direction (solid arrows in figure I7.9P . H 

Then, for large enough forces, particles are more energetic and there is more proba- 
bility of being launched over the barrier by thermal fiuctuations; then the wave-refiection 
dominates and we see an amplification of the rectification by quantum effects. On the 



^ Recall that in quantum mechanics, much as particles with energy lower than the barrier can be 
transmitted to the other side, by tunneling, the complementary effect, equally counter-intuitive, can also 
take place: particles with energies larger than the barrier can bounce back: overbarrier wave-reflection. 



classical 




85 



CHAPTER 7. APPLICATIONS I: BROWNIAN PARTICLE IN PERIODIC POTENTIAL 






Figure 7.9: Saw-tooth potential in the absence of force, and as deformed by zizF. The 
barrier is fower to the easy side than to the hard side. The soHd arrows indicate the 
transmission for energies smaller than the barrier (tunnelling) and the broken arrows the 
transmission for energies larger than the barrier (affected by overbarrier wave reflection). 
The length of the arrows indicates qualitatively the probability of those processes. 

other hand, for small forces/energies, it is more probable to find particles with energies 
lower than the barrier. Here tunnel events win, and the rectified current is reduced. 

In PAPER I we presented results for fixed amplitude F, studying the dependence of 
the current on the temperature. 

7.5 example III: AC response in periodic potentials (LRT) 

In the previous examples we have addressed stationary transport in a cosine potential, 
subjected to a constant force, and then rectified currents in a ratchet potential, with a 
driving force in the adiabatic regime. Thus, we could reduce the studies to the computation 
of stationary solutions, and combinations of them. Now we are going to address a genuine 
AC problem, with the driving out of the adiabatic regime. 

Let us return to the cosine potential (j7.9p . and let f{t) = 6- e"^*. If the amplitude is 
small 5 <ti 1, we can use linear response as in section [HTTl The observable of interest now 



To compute the amplitude x('^) we can use the continued-fraction method as explained 
in section 16.5.21 Note that, being in linear response, {p) is proportional to the force, so 
that {f{t)) = gives a zero time average of the current, irrespective the potential. Then, 
the asymmetry of the potential does not bring much new, and for simplicity we return to 
the cosine potential. 

Figure 17.101 shows details of the response in the low frequency range (left panel) and 
at high frequencies (right). The low frequency spectrum shows a Debye profile, while at 
high frequencies we see a pattern of absorption lines (we just plotted the imaginary parts). 
These figures can be understood qualitatively recalling the discussion of section [Ol There 



is the AC current [cf. Eq. (i630|) ] 




(7.20) 
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Figure 7.10: Left: low frequency range of the susceptibility x(w) (cf. figure [6^21) . The real 
and imaginary parts are plotted for K = 1000 (classical), T = 0.5, 7 = 0.1 and F = 0.1. 
Right: x{'-^) ill the high frequency range (imaginary part only), for K = 200, T = 0.05, 
and F = Q (as in PAPER I). The 7 = 0.01 curve is essentially the classical result, and 
for 7 = 0.0003 and 7 = 0.0001 the peaks become narrower. The vertical lines mark the 
transition frequencies of the associated approximate quartic oscillator (see the text). 




we argued that in general the response would be a combination of Debye's and absorption 
curves, pretty much as what we see in figure I7.10[ The absorption peaks are Lorentzian 
curves centered on the energy differences, and with a width oc 7 (or inversely proportional 
to the decoherence times, if you prefer). On the other hand, the broad relaxation maximum 
of the low frequency Debye curve (note the log scale), would give the relaxation time. 

The reason of the Debye part is the following. Along with the erasing of the non- 
diagonal terms in the density matrix (decoherence), the system has to readjust the popu- 
lations of the wells (the cosine potential is now tilted by the AC force) . This leads to the 
relaxation time (slow dynamics) we find in the low frequency range. 

As for the high frequency curves, the widths of the absorption peaks can be cal- 
culated semiclassically (PAPER I, appendix F). We performed a WKB calculation of 
the widths of the lower bands for the example displayed {K = 200), finding = 
(4/if) exp[— i^(l — e)/2] with e = E/Vq. Therefore the widths are exponentially small, 
and for K = 200 they can be considered as discrete levels (recall Fig. 12. 3p . Besides, their 
position can be estimated from the energy levels of the anharmonic oscillator obtained 
expanding the cosine to fourth order cos(x) = ax^ + bx^. Indeed, plotting the associ- 
ated energy differences, we find that the absorption peaks are centered around them, as 
expected. Now, if we increase the damping, we see how the absorption curves broaden, 
with a width proportional to 7 oc I/T2. Eventually, as the damping increasingly blurs the 
energy levels, the curves approach a smooth peakless "classical" profile. 



87 



CHAPTER 7. APPLICATIONS I: BROWNIAN PARTICLE IN PERIODIC POTENTIAL 



7.6 summary 

The example of the particle in a periodic potential has allowed us to review many of 
the concepts and tools introduced in the preceding chapters: Fokker-Planck equations, 
quantum master equations, phase-space methods, relaxation and decoherence times, . . . 

The continued-fraction technique, plus phase space, has made possible to solve the 
problem as an extension of Risken approach to the classical problem. In this way we 
have bypassed the problem of the underlying spectrum (demanding, as is continuous and 
non-bounded). During the whole chapter we had the classical limit as reference, in order 
to understand our results as corrections (recall that the Caldeira-Leggett equation is 
semiclassical). 

We finally remark the differences between the classical and quantum distributions 
(recall the wave-mechanical interpretation). We have discussed how the interplay of over- 
barrier wave reflection and tunneling can increase or decrease the rectified current in an 
asymmetric potential, depending on the force/energy range. We have also seen how the 
damping can make the curves to approach the classical counterparts. We closed the chap- 
ter with the study of the dynamical susceptibility spectra and characteristic times. 
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Chapter 8 



applications II: 
statics &; dynamics of 
superparamagnets 

Quantum nanomagnets are an example of mesoscopic objects described effectively by 
"one" degree of freedom (their net spin). In this chapter we are going to study the statics 
and dynamics of these systems, where the interaction with the environment is important. 
We start reintroducing the model and investigating its equilibrium properties. Next we 
specify the master equation to be solved to study non-equilibrium behavior. As examples 
we will address the relaxation mechanisms (analogous to the calculation of Ti with Bloch's 
equations), by analyzing the longitudinal susceptibilities (linear and non- linear). We will 
conclude with the calculation of the transverse susceptibility (study of T2). 

8.1 spin Hamiltonian 

Superparamagnets are solids or clusters of nanometric size where each molecule can be 
described by a spin Hamiltonian |21j : 

= -DSl - B,S, . (8.1) 

Here D is the anisotropy parameter, originating from the spin-orbit coupling, while —BzSz 
is the Zeeman term accounting for the coupling to external fields. As a first approximation 
the magnetic molecules can be considered as independent of one another, so the total 
Hamiltonian is a sum of contributions like (j8.ip . Its eigenvalues are e„ = —Dm? — Bzm. 
For D > the anisotropy term in ()8.ip yields a bistable structure with minima at 
m = its, while the Zeeman term tilts this "double well" (figure 12.41 chapter [2|). In 
equilibrium the spin orientations are like those of a paramagnet (the magnetization and 
susceptibility are given by Brillouin and Curie type laws), but due to the large S they are 
called superparamagnets. Physical realizations of this model are the Mni2 or Feg clusters, 

^ In this chapter we set h — fiB = ks ~ 1; remember the footnote [S] in chapter [5] 
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among others [22], which have a moderately large spin S ~ 10, and where quantum effects 
of the spectrum's discreteness can be studied [1]. 

It is convenient to introduce scaled variables (/3 = 1/T): 

a = (3DS\ ^ = PB,S, h= \ . (8.2) 

The parameter h is called the reduced field. When there is no external field obviously 
h = 0, whereas h = 1 just when the barrier disappears. In these variables, the thermal 
Hamiltonian reads PTis = —a{m/S)'^ — ^{m/S) and we will take the classical limit 5 — > 00 
keeping a and ^ fixed (note \m/S\ < 1). In this way the energies are kept constant, while 
introducing more and more levels towards a continuum. Thus the spin number S plays 
the same role as K in the particle problems (characteristic action over h). In that case 
when K ^ 00 the number of bands below the barrier was increased, but the potential 
height was kept fixed. 

Prom a formal point of view, varying parameters allows the study of problems with 
equispaced spectrum {D = 0), or problems with more featured spectra {D ^ 0). In the 
case > 0, we have mentioned that the spectrum has a double-well structure, while 
D < has a single well (figure [27i]) . The level spacings are 

em - em±i = Amm±l = ± [^(2m ±1) + B] (8.3) 
so that, when Bz = 0, the levels m and —m become degenerate. 



8.2 Hamiltonian coupling to the environment 

As these spins "live" in a crystal, they interact with the lattice phonons, providing an ex- 
ample where the coupling to a bath of oscillators is well justified. Thus we write [cf. (|3.26|) ] : 

Htot =ns + Y^ Fa{S) (g) c„(6+ + b.a) + ^c^b+ba (8.4) 

a a 

where TCs is the Hamiltonian accounting for renormalizations (recall sect. 13. 4| ). As we see, 
these systems can be studied within the open-system formalism of the previous chapters. 
The coupling between the spin and the phonons includes 

FiS) = r]+[v{Sz), S+]^ + 7^4v{Sz), 5_]^ , (8.5) 

where the complex conjugate constants ?]_ = r/!j_ ensure hermiticity, the anticommutator is 
represented as [ , ]+, and v{Sz) is any function of the operator Sz- Note that the coupling 
is realized through the ladder operators S± = Sx^iSy, which obey S±\m) = im,m±i\'m±l) 
with the custom factors im,m±i = \/S{S + 1) — m{m ± 1). 

In the case v{Sz) = cte, the coupling F oc r]^S^+r]^S+ is linear in the system variables. 
This is analogous to the coupling used in the particle chapter, and here will serve us as 
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control model. On the other hand, v{Sz) oc Sz for spin-phonon (magneto-elastic) coupling 
[85\ [86], and we have to use the more realistic F ~ [5*2, S"-!-]-!-. In the former, linear 
coupling case we use an Ohmic distribution of bath oscillators, J{u!) = \lo, as we did 
previously. But when studying the magneto-elastic model we will use a super-Ohmic bath 
with a spectral distribution J{io) = Xco^ appropriate for phonons (or photons) in three 
dimensions [I]. 

Both in the equilibrium properties and in the dynamics we will need the explicit form 
for the matrix elements Fnm = {n\F\m): 

= r/+[w(m) -I- v{m!)Yrn,m' ■ (8.6) 

Note that when v{m) ^ cte (spin-phonon) the L.^^^' add an extra dependence on m in 
the coupling/damping, compared with the linear case; if you wish, the spin analogue of 
position-dependent damping in particle problems |14j . 



8.3 equilibrium properties 

In chapter H] we discussed the equilibrium properties of open systems. In the previous 
chapter we already showed stationary Wigner distributions for a particle in a periodic co- 
sine potential. There, as we worked with 7/T <C 1, the stationary solution was essentially 
the canonical distribution at zero damping (figure 17. 3p . Now we want to study explicitly 
corrections to the thermodynamical properties when A 7^ in spin problems. 

Let us recall that to obtain the damping-dependent corrections we had to compute 
certain Z^'^h Then Ztot/^h — -2^s + -2^^^^ with ^toti and Zg the partition functions of 
the whole system, the bath and the closed system, respectively, Eq. (j4.9p . From Z we can 
obtain the free energy, and hence all other thermodynamical quantities, as in section! 



For a spin-bath coupling linear in S± as in ()8.6p . we have |FnmP = ^n+i n^mn+i + 
(^-n n-i^m,n-i- Plugging this l-FnmP in expression (j4.10p for Z^'^'> we have 



2 



e^ 

m 



^m+l,m^ (Am+l,m) + ^m,m-l^ {^in-l,r 



(8.7) 



where we have used the energy levels em = —Dm?—Bzm, while Am,m±i = i [D{2m ±1) + Bz 
are their differences ("transition" frequencies). 



8.4 example 0: equilibrium of the isotropic spin 



If we consider D = (isotropic spin) the spectrum becomes equispaced and the arguments 
of the <I> in ()8.7|) are simply ±Bz (the constant level difference). As these do not depend 
on the index m, they can be taken out of the sum, as in the example of the harmonic 
oscillator in section Then Z^'^^ can be explicitly calculated [cf. Eq. (j4.2i 



2T 



+ coth 



^^ = ^{Bz)±^[-Bz 
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Here Z^^ is the partition function of an isotropic spin Z^ = sinh[(5'+l/2)i?^/T]/ sinh[i?2/2T] 
and Is/lf^ the corre 
Brillouin function: 



and M^^^ the corresponding magnetization (at zero couphng), which is nothing but the 



Mi°) = (5 + i)coth[(5 + i) 



1\ 

T 



1 ^ 
- cotli 



IB, 
2 T 



(8.9) 



Now, with our Z^'^^ we can compute the magnetization including the corrections due to 
the bath couphng: 

d Z(^) 

^ Mf) + T-—-— . (8.10) 
(JBz Zg ^ 

Now we take the derivative and plot the results in figure ISTTl These are the corrections to 
Brillouin's law for several spin values, and as a function of the field and the temperature. 
On the left panels we plot the bare Brillouin curves as reference. 

But let us first discuss the isolated spin limit. At zero temperature the system is in 
the ground state m = S and hence the magnetization is "saturated" = S. Increasing 
T the levels S — 1, 5 — 2, ... become populated, and starts to decrease. Eventually, 
at very high temperatures all levels are equally populated, and the magnetization goes 
to zero, say, by thermal misalignment from the field direction [M, = Yl'^^Qmrm and the 
T ^ oo population of the levels m and —m is the same). 

Keeping these A — > results in mind, let us turn to the corrections due to the bath. 
First, we see that the finite-coupling corrections vanish in the S ^ oo classical limit, and 
in the high temperature limit too, as we could have expected (section [4.2p . Besides, the 
corrections are negative, and they are more noticeable at intermediate values of B^/T. H 

(2) 

It is worth remarking that at T ^ the correction remains finite, Mr' < 0, and 
hence = M,^^ + M^^^ < S, contrasting the saturated Mz{X = 0) = S*. In other words, 
switching on A 7^ the excited levels became populated. What should we think about 
this? Well, we have to take into account that the system proper is the total one, spin 
plus bath. Besides, [Hs + 'Hh,'Hsh] / 0. Then the ground state is not \S) |vacuum) 
(with I vacuum) the bath's ground state), but spin and bath are entangled. Seen in other 
way, the reduced density matrix is not diagonal in the basis, which makes Qmm foi' 
excited states even at T = 0. 

We already mentioned in chapter H] that in the systems we have in mind this type of 
corrections are probably small, and we do not expect they being very relevant in a Mni2 
or Fes, for example. Maybe very low temperature measurements, to check < S. This 
would give a quantitative measure of the entanglement between spin and bath (see Jordan 
and Buttikker \87\ [88] for a similar discussion in the harmonic oscillator) . Anyway, though 
the equilibrium formalism helped us to understand and correct some misconceptions, in 
what follows (dynamics), we will consider for all practical purposes that the equilibrium 
state is the closed canonical one. 



We obtained results qualitatively similar for the anisotropic spin with D > 0. 
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Figure 8.1: Top panels: Brillouin function (left) and its dissipative corrections (right) for 
several S as a function of = B^SjT . For the corrections we use AS* = 0.05 (corresponding 
to the classical damping, see footnote [3] below) and J{oS) = Xlo, in the model linear 
in S (section 18. '2p . Bottom panels: the same quantities (left Brillouin, right dissipative 
corrections), but plotted vs the temperature (complementing the top plots vs. 1/T). 
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8.5 dynamics 

It will not surprise the reader if to study the non-equilibrium properties we use the quan- 
tum master equations of chapter O Specifically, we will use the second order master 
equation (|3.4ip . in what we called the improved rotating- wave approximation (RWA) 

jim/dt — i^nmQnm ~l~ T^nn+l,mm+lQn+lm+l ~\~ T^nn.mmQnm ~\~ T^nn—l,mm—lQn~~lra~l ■ 

(8.11) 

At the end of this chapter we will briefly touch the comparison between crude RWA, and 
its improved version. 

At some point we would like to add a small transverse field, e.g., B^^ to study the 
transverse response Xx{^)- Then the system Hamiltonian 7is becomes 

= -DSl - B,S, - i(S+5_ + (8.12) 

with B± = Bx HBy. But for arbitrary B± we do not know the eigenvalues and eigenvec- 
tors of Tis analytically, so we cannot calculate the unperturbed evolution of the coupling 
term FnmiT) = -^rime"'^"™^ (which was required to calculate the relaxation term TZ in 
section [3.5. ip . However, restricting ourselves to B± ^ 1 we can approximate the required 
evolution by the diagonal part. Naturally, we can and will fully keep the transverse field 
in the unitary evolution. All things considered, we will have to solve 

d^'nm/dt — l^nmQnm 

"I" ~\~ ^mm—lQn—lm) ~\~ ~\~ ^mm+lQn—lm) 

~l~ 'TZnn+l,mm+lQn+lm+l ~^ 'TZnn,mmQnm ~\~ 'T^nn—l,mrn—lQn—lm—l 1 (8.13) 

where A„m are the energy differences for the longitudinal part (18. 3p . For the explicit form 
of the 7^'s, see PAPER IV. 

A final comment on the resolution method. Note that the master equation ()8.13p . 
can be formally written as Qnm = Yla i3 Q^^+a Qn+a^m+p- Then, introducing the matrices 
[Qn,n+a]m,m+/3 = Qn'n+I and the vectors [C„]m = Qnm, we have 

Cn = Qn,n-lC„_i + Qn,nC„ + Qn,n+lC„+i . (8.14) 

That is, we have rewritten (j8.13p in a form suitable for the use of continued fractions 
(chapter El PAPER IV). 

8.6 examplel: longitudinal relaxation &; linear susceptibility 

Our first example of dynamics is a study of the relaxation mechanisms by means of the 
linear susceptibility. We will consider that no transverse fields are applied {B± = 0), and 
analyze the time evolution of = {Sz): 

S 

Mz{t) =Tv{SzQ) = Nm = Qmm- (8.15) 

m=-S 
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In this problem we just need the diagonal density-matrix elements to construct the re- 
sponse Mz- Besides, for B± = the diagonal and off-diagonal elements in (|8.13p get 
decoupled. Then it suffices to solve the system of 25 -|- 1 coupled differential equations for 
the populations Nm {balance equations): 

= {Pm\m+l^m+l — Pm+l\m^m) + {Pm\m-l^m~l — Pm-l\m^m) (8.16) 

where the transition probability is Pm\m' = T^mm' ram'- I 

We will calculate the linear response of to an ac perturbation 6Bze~^^^. To this 
end, recall that Xz{i-^) could be constructed as in Eq. (|6.3ip . using the eigenvalues and 
eigenvectors of the balance equation without perturbation; and that equation (16. lip relates 
the response in the time and frequency domains. 

We can compute Xz{^) directly with continued fractions too, as in section [6.5.21 How- 
ever, the eigenroute not only gives Xzi^), but allows to analyze and characterize its con- 
tributions (while continued fractions are a kind of a black box). 

8.6.1 eigenvalues and eigenvectors of the relaxation matrix 

Let us first write the balance equation (|8.16p in a matrix form, with the vector of popu- 
lations [N]m = Nm 

dN/d^ = -7^N. (8.17) 

The eigenvalues Aj and eigenvectors |Aj) of the {2S -|- 1) x (2S + 1) matrix TZ can be 
obtained by numerical diagonalization. We have plotted them in figure 18.21 One of the 
25-1-1 eigenvalues is zero, corresponding to the stationary solution. Then we can write 
the solution of (j8.17p as N = Yl'i=i e"^** |Aj) + |Ao) with | Ao)m oc exp(-/?em), that is, 
the canonical distribution (second panel of figure 18. 2p . 

The remainder eigenvalues are positive (ensuring convergence to the stationary solution 
as t ^ oo) and real (the density matrix is Hermitian). Figure [821 also shows that the 
eigenvalues are organized in two "groups". Far away from the others we find Ai, which 
generates the slow long-time dynamics (ri = A^^^). Indeed the structure of the eigenmode 
|Ai) (figure 18.21 3rd panel), with population decrease in one side and increase in the 
other, corresponds to transfer from one well to the other. As tunnel is not present (no 
transverse field), this population transfer should be channelled over the barrier top, and 
hence we call it the overbarrier process. On the other hand, the rest of eigenvalues give 
population readjustments inside each well (fast dynamics); we call these the intrawell 
processes. The former, Ai, depends exponentially on the barrier height, whereas the 
intrawell modes depend only polynomially. Whence the large separation of the scales of 
these two processes seen in figure 18.21 (note the logarithmic scale) . 

^ As we work here with the diagonal elements only, while Qmm is the probability of having the spin in 
the state m, we can take the classical limit directly. On so doing, we get the longitudinal Fokker-Planck 
equation (I2.23p . In the general case (with non-diagonal terms), we have to resort to phase-space methods, 
as in chapter [5] Taking the classical limit allows us to relate both theories, classical and quantal. In 
particular, for Ohmic baths, we get the relation AS = All where A is the damping used here, whereas All 
is the phenomenological damping in the Landau-Lifshitz classical equation (see PAPER III). 
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Figure 8.2: Eigenvalues and eigenvectors of the matrix associated to the system (I8.16p . 
The parameters are a = 15, = and S" = 10 . The couphng is magneto-elastic and 
J(a;) = Xu!^, with A = 10~^. Top panel: non-zero eigenvalues. Lower panels: a few 
representative eigenvectors: |Ao) (stationary solution); |Ai), determining the overbarrier 
dynamics, and two fast modes IA17) & jAig), associated to the dynamics in the right and 
left wells, respectively. 



96 



8.6. EXAMPLE I: LONGITUDINAL LINEAR SUSCEPTIBILITY 



Well, having the eigenvalues and eigenvectors p we can construct now the dynamical 
susceptibility. The eigenvalues, being real, would produce a susceptibility in the form of 
a sum of Debye profiles (section 16. 3p 

We have plotted a few susceptibility curves in figure [831 The response is given, at most, 
by two main contributions. We associate the low frequency peak with the slow overbarrier 
process. While the peak at high frequency, when present, is to be attributed to the fast 
intrawell dynamics. 



8.6.2 bimodal approximation (analytical) 

To analyze the susceptibility curves, given the mode's structure above, we tried a relaxation 
profile consisting of two effective modes (as Kalmykov et al. did in the classical case [89]) 



eq 



ai 1 — Ol 

+ 



1 + i UJTi 1 + i IjJT^ 



^.18) 



Here x^'^ is the equilibrium susceptibility and ri = A^^ ^, which follows an Arrhenius- 
Kramers law: 

ri oc e^^^ (8.19) 

This exponential dependence on the barrier height A{7 is typical of problems of escape 
out of a metastable minimum [901 191j . The specific calculation for our problem is given in 
PAPER III. 

The susceptibility above has two more parameters, ai and r^, which can be determined 
following the idea proposed by Kalmykov et al. in the classical problem [89]. They argued 
that we have two parameters to determine, while Xz{'^) can be written in closed form in 
both the limits of low and high frequencies; in terms of Tint for w — > [Tint is the area 
below the relaxation curve, Eq. (|6.32p ] and in terms of Tpf at high to [the initial slope of the 
decay ()6.33p ]. The advantage is that both Tint and can be computed analytically. Then, 
these two limits would provide the two unknown parameters oi and t^, so yielding a closed 
formula without fitting parameters (PAPER III) . The results of the analytical ansatz are 
the lines in figure 18.31 which account quite well for the symbols, obtained constructing 
Xz{'^) from the numerical eigenstuff. 

Let us comment on the curves. The relative importance of the two effective modes 
is governed by the parameter oi: when oi = 1 there is only overbarrier response, while 
switching to ai = 0, the whole weight is shifted to the intrawell response. Figure [831 shows 
how this is realized by changing the field Bz, with ai = 1 at zero field and ai decreasing 
towards zero at large fields. Indeed, at large h there is no barrier, so the overbarrier 
contribution is naturally expected to vanish. However, figure 18.41 shows that ai = much 
before the barrier has disappeared. 



* To avoid possible misunderstandings, let us make clear that these eigenvectors and not vectors of the 
Hubert space of the problem, but they are rather n-tuples containing the diagonal elements of q. 
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Figure 8.3: Real and imaginary parts of Xzl"^) for a spin 5 = 10 in various fields h. The 
parameters are a = 15, T = 0.1 and A = 10^^ (magneto-elastic model with J{ijj) = Xuj'^). 
Note that the frequency is scaled by ri. The symbols are the numerically exact results 
and the lines the bimodal approximation (18.180 . 
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Figure 8.4: Coefficients (k = 1, 2) plotted vs. h. The ak are approximately independent 
of tlie spin-batli coupling considered (they depend on tlie parameters of tfie spin Hamil- 
tonian and T). Here we used 5 = 10 at o" = 14. Notice that while the barrier disappears 
for h = 1 (Sec. IS.ip we have = quite before that value. 

Indeed, figure 18.31 shows that at /i = there was only one peak, associated to the 
overbarrier mode. Physically, at zero field the fast dynamics of the two wells is equivalent 
(see figure [8^ lowest panel), and it averages to zero when looking at (Sz)- At intermediate 
fields h = 0.15 the response clearly exhibits two modes. Increasing h further only the 
intrawell peak is left (though there is still barrier at /i = 0.3, the upper well is thermally 
depopulated). 

We close with the isotropic spin. If we set D = , we always have a single peak (see 
[92j and PAPER IV). In a sense, this is clear now for us; in the isotropic spin there is no 
barrier imparting a large time scale separation, and this kind of relaxation measurement 
does not resolve the differences between the fast modes. 

8.7 example II: non-linear susceptibility 

Let us illustrate, with the example of the longitudinal relaxation, how to proceed if we 
want to go beyond the linear response regime. In this and the next section we will see how 
going to higher orders in the response, one can obtain more information on the system 
|93l I94j . In the nonlinear case, we do not have anymore an equation like (16. lip relating 
the behaviour in the time and frequency domains. Here we are going to study directly the 
response as a function of the frequency of the driving field. 

To get the dynamical susceptibilities, we excite our system with 5Bze^^^ and proceed 
as we did in section [6.5.21 That is, we start from the master equation written formally as: 



and we Fourier expand the evolution operator C = Cq + 5 ■ e'^*vCi + (5^ • e'^'^*£2 + • • • as 
well as the solution Q = Qq + 5 ■ e^^^Qi + 5^ • e^^'^*£>2 + • • • Here we have explicitly extracted 



dtQ = Cq 



(8.20) 
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the temporal dependencies e'^^'^* from q and C. Therefore, Qk and Ck {k = 1,2, ... ) are 
time independent (of course, qq and Cq too). Now, equating order by order we have 

= CoQo (8.21) 
iiogi = CqQi+CiQq (8.22) 
2\ujQ2 = C0Q2+ C\Qi + C2Q0 (8.23) 



(we have written to second order only, but the perturbative structure is clear). Once 
again, this chain of equations can be solved by continued fractions, as in 16.5.21 First we 
solve for qq; this is then inserted in the second equation, and we solve for qi; knowing qq 
and Qi we can solve for Q2 in the third one. 

In some cases this set can be solved analytically. For example, for spin ^ we had the 
equation of motion for the magnetization given by Bloch equations (section I6.2p 



= -TrM, + M^{oo) (8.24) 

Then £ = F^ (a number) and operating like in the chain ()8.2ip - ()8.23p we find simple 
expressions. To first order we recover the Debye (|6.24p . while for the first nonlinear 
susceptibility we have 

/ A eq 1 eq ^^t[ , _ Bti 

X2(w = X2V— ^1 TTT^ vn~r~^^ ^ ' '^1 = • ^'^^ 

l + \lujTl [1 + \UJTi){l + \ luJTi) OBz 

Here we have written ri for the relaxation time, while = d^Mz/dB^, k = 1,2 are the 
equilibrium linear and nonlinear susceptibilities. As we see, the nonlinear response gives 
information, not only of ri, but of its iJ^-derivative as well, t[. As mentioned above X2i^^) 
seems more sensitive to some parameters of the problem. 

Well, this was for S" = ^, but we want to get the response for any S. In general, the 
response will include relaxation "blocks" of the type (I8.25p . and maybe cross-terms. As we 
did for the linear susceptibility, we will try to find an approximate analytical modelization. 



8.7.1 generalization of the bimodal approximation to nonlinear response 

The reader could rightly guess that we are going to try a similar reduction of the 25" modes 
to 2 effective ones (intrawell and overbarrier) to account for the nonlinear response. Thus, 
we construct an effective model with dynamics governed by the two relevant relaxation 
times of our problem: ri and Tw- Although now we need to use 2x2 matrices in (|8.2ip - 
()8.23p . one readily gets: 

eq «2 eq aiiujT[ 



X2{UJ) = X2 . , -o Xi 



l + i2wTi (1 + ia;Ti)(l + i2wTi) 

eq (I-Q2) eq (l-aQic^T; 

l + i2wTw {l + iuJT^){l+\2uJT^) ^ ■ ' 
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Figure 8.5: Real part (squares) and imaginary part (circles) of the linear (top) and 
nonlinear susceptibility (bottom). The symbols are the numerically exact results and 
the lines the bimodal approximations (18.180 and (18.250 . We have set S = 10, and the 
parameters a = 14 and ^ = 5 in reduced units ()8.2p . The coupling is magneto-elastic with 
F ~ S,S±. 




where 

t[ = dn/dB, r; = dr^/dB, (8.27) 

This is a generalized bimodal formula for the nonlinear response, without any free param- 
eter, as before. The new coefficient, a2, tells when the intrawell modes enter the scene 
(decreasing from 1), and was already plotted in figure EH We see that for a given field 
a2 > oi, so that the fast modes cannot be detected at lower fields using X2i^)- Had we 
obtained a2 < ai, we would have used X2('^) to detect the intrawell modes without needing 
as large fields as when using xi- 



101 



CHAPTER 8. APPLICATIONS II: SUPERPARAMAGNETS 



In figure 18.51 we compare the new formulas with the numerical results obtained by 
solving the chain (j8.2ip - (|8.23p with continued fractions. We see how the bimodal ansatz 
works very well in the nonlinear case too. We also see how the intrawell modes as less 
noticeable in the nonlinear susceptibility X2i^), as we anticipated from 02 (-Bz). 



8.8 example III: ri and tunneling — experiments in Mni2 

In this section we will use the nonlinear formalism for the analysis of experiments in the 
superparamagnetic molecular cluster Mni2. In particular, we will see how to get more 
information on the relaxation time ri from the X2(<^) vs. lo spectra. 

Well, in actual superparamagnets the Hamiltonian is a bit more complicated than (jS.ip 

Hs = -DSl - B,S, - E{Sl - SD - B^S, + • • • (8.28) 

The term —E{S'^ — Sy) is a consequence of the distortion of the anisotropy axis, while 
Bx accounts for possible transverse fields. But when these terms are small they can be 
treated perturbatively. 

Let us first discuss (with the hands) how these extra terms, which do not commute 
with Sz, affect the relaxation time ti (I8.19p . In the absence of transverse terms, the states 
\m) and \—m) are degenerated at Bz = 0. Now, the terms E and Bx mix \m) and \—m) 
(figure 18. 6p . The new eigenstates can be formed from their symmetric and antisymmetric 
combinations |s) and \a) 

\s) ^ -^(|m) + |-m)) 

\a) - -L(|„^)-|-rr^)) (8.29) 

with e\a) ~ £|s) = ^ the tunnel splitting, which can be calculated perturbatively p6l [95] . 
These states are delocalized (in m), and the spin can "tunnel". 

If we now apply a field Bz the degeneration is lifted and we block the tunnel channel. 
This can be seen from the perturbed eigenstates: 

1/2 

Q+|m) + a_| — m) 

The transverse terms are assumed small ($7 ^ 1), so we have a+ = 1 and a_ = at not 
very large fields Bz- Therefore, out of resonance we have \s) = \m) and |a) = | — m), and 
we are back in the behavior with no transverse terms. 

However, increasing the field up to Bz = D, the levels \—m) and |m — 1) enter in 
resonance, making possible tunneling through these states. This is repeated at all crossing 
fields, Bz = kD, where the states \—m) and \m — k) become degenerate (figure [2^ . 
Therefore, by changing the field Bz, we "open" and "close" tunnel channels. 

The tunnel splitting decreases exponentially from the states near the barrier top to 
the lower laying states near the bottom of the wells. This can be understood because the 



a± 

a-\m) — a+\ — m) 



(8.30) 
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Figure 8.6: Schematic representation of the sphtting of the levels |m) and \—m) under 
resonant (top) and non-resonant conditions (bottom) as discussed in the text. 



terms S"^ and Sx couple states with m it 2 and m it 1 in each perturbative order: then, to 
connect S all the way to —S we need to go to higher orders in the perturbation. H 

On the other hand, in the temperature ranges we are going to address (T > 2 K) tunnel 
proceeds incoherently. That is, the decoherence time is so short, that coherent oscillations, 
left-right, are never observed [93]. In this incoherent regime, the relaxation time ti can 
be written as [1]: 



'1 



"' = '^r^,n^'e-^(^-'"-^-^) (8.31) 



where the sum extends to the resonant pairs m and m' . As we said, for the lowest laying 
pairs Tm^m' — 0. As a matter of fact, the barrier that the spin has to overcome decreases 
in the crossing fields (as the tunnel channels get open) and we have ri oc e^'^^'=f with 
AC/gf < AC/. In figure 18.71 we have plotted the formula (j8.3ip for a simple model with 
transverse field and displayed experimental data for Mni2. At the crossing fields, where 
the levels enter into resonance, the relaxation time shows dips, as a consequence of the 
decrease of the effective barrier |95] . 

How could we see this behavior studying X2{^) curves? Well, if we work with not too 
strong static fields, we have ai = a2 = 1. Besides, if we probe the system with frequencies 
of the order of r^"^, the fast intrawell contribution to ()8.26p is negligible {ujt^ <^ 1, since 
T„ ^ Ti). Then we can work with an approximate X2('^) dominated by the overbarrier 
contribution 

X2(W) = X2 \ , _ -Xl 



l + \2ujTi (l + iLjri)(l+i2a;ri) 

At fields Bz close to the crossing fields, the derivative r( becomes large (we are inside 
the dips). Therefore, we expect that the term with r( would dominate in X2('^) in those 
ranges. On the other hand, the derivative changes sign, and so would do x^il^) at the 
crossing fields. 



® This agrees with the semiclassical picture where tunnel is smaller the higher the barrier over the 
tunnel channel. 
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Figure 8.7: Left: the relaxation time ri from Eq. (j8.3ip . The curve corresponds to 
including a transverse field in "Hg = —DSl — BzSz — B^Sx with = 10~^, at a = 15 for 
a spin S = 10. We have also plotted for reference ri for Bx = (thin line), where there is 
no tunneling, and the relaxation time does not display dips at the resonance fields. Right: 
Ti in Mni2 at T = 8K from susceptibility measurements (PAPER V) . 
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Figure 8.8: Left: prediction for X2i'^) (imaginary part) with the theoretical ti from 
figure 18.71 left, plugged into the analytical (|8.25|) . Right: experimental X2{^) curves 
corresponding to the ri shown in figure \87l\ right. 
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We have plotted X2i^) as a function of Bz in figure EHl The experimental data shows 
the same qualitative features as those just discussed on the basis of the formula (j8.32p . We 
thus see that X2i^) is a good tool to detect tunnel in these systems, as it provides resonant 
features (peaks and sign changes) not easily mistaken. And all this because X2{^) is more 
sensitive to the variations in ri that the linear susceptibility (PAPER V) . 



8.9 example IV: transverse response 

The linear susceptibility is still useful, as it highlights some other aspects of the problem. 
In this last example we are going to compute the transverse susceptibility, as we already 
did for the spin ^ case. To this end, we excite with a transverse perturbation 6 • B^g 
and we monitor the response of = {Sx) 

= 2 ^'^+ ~'~ "^"^ ~ 2 {^m-l,mSm-l,m + im+l,mQm+l,m) (8.32) 
m 

In this problem, the unitary part in the master equation ()8.13p mixes off-diagonal and 
diagonal elements, incorporating the ladder factors £m,m±i = \/S{S + 1) — m{m ± 1). 
Therefore, we have to consider all the density-matrix elements. This can be done solving 
with the continued fraction method (section l6.5.2p . The response obtained is displayed in 
figure El 

Based on the results of the spin-1/2 in 16.2.21 one could expect some absorption spec- 
trum. The imaginary part is indeed made of Lorentzians centered at the level differences. 
But due to the anisotropy, the Tis spectrum is not equispaced now, and we can find up to 
25 + 1 peaks (cf. the cosine problem in figure [7TTU|) . 

Following our custom, we will try to draw an approximate formula over the exact 
curves. Recall that under the crude rotating- wave approximation Eq. (j3.40p . the time 
evolution of the non-diagonal elements is simply given by Qm±i,m oa exp[— i(Am_,_im + 
T^m±im±i,mm)t]- Then the response in the frequency domain follows from Eq. (j6.1ip as 



'i(u; + Am-lm) + T^m-l,m i{iO + Am+lm) + 



.33) 



Here Xx^ is the equilibrium susceptibility, while TZm±i,m = T^m±im±i,mm are relaxation 
coefficients and the amplitudes depend essentially on the population difference between 
the states m and m + 1. 



_|_ 1 Q~P^m g~/3em=Fl ^ 

'^m 9 eq ^ T ^m.m^l • (8.34) 



This constitutes a generalization of equation (j6.28p for S > ^. This formula is plotted as 
the continuous line in figure 18.91 where we see that it fits quite well the numerics (which 
correspond to the improved RWA, Eq. (|8.13p . symbols). 

If we now increase S or the coupling A (that is, decreasing the quotient Amm±i/^), 
we have observed that the crude RWA (with which we got the formula) starts to fail 
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Figure 8.9: Imaginary part of the transverse susceptibility Xxi^) for a spin S = 10, at 
a = 5 with a weak damping A = 3 • 10~* (magneto-elastic couphng). Only the positive 
frequency is shown, as Xxi^) — ~Xxi~^)- The symbols are the exact continued- fraction 
results (PAPER IV), and the continuous line is Eq. (j8.33|) . We have also marked the level 
differences with vertical dashed lines. 

(see also ^37j for more RWA failures). We could still try to get the classical limit of the 
susceptibility [96]. But there the formula would be valid only in the limit A = 0, anyhow 
recovering Gekht's formula for the absorption line of a classical superparamagnet |97i 
I96j . In any case, for the couplings and spin values in ordinary quantum superparamagnets, 
equation (j8.33p could be a valuable tool to fit experimental absorption curves. 

8.10 summary 

In this last chapter, we have sought to apply all the formalism discussed along the thesis to 
the example of quantum superparamagnets. What we learnt in chapter U] has been used to 
calculate dissipative corrections to the celebrated Brillouin law of quantum paramagnets. 
With the master equation of chapter [3] together with the methods presented in chapter El 
we have studied relaxation mechanisms in superparamagnets. We could contribute to 
the study of the phenomenon of tunnel of the magnetic moment, studying the nonlinear 
susceptibility. Finally, we derived approximate expressions for several quantities; the 
formula for the dynamical transverse susceptibility can be used to estimate the decoherence 
times in superparamagnets. 
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As it is the custom, let us conclude remarking the most important results presented in 
this thesis. 

• In the classical limit we worked with the bath-of-oscillator Hamiltonian ()2.4p to 
describe open systems, obtaining the corresponding Langevin and Fokker-Planck 
equations. The model allows to understand the paradoxes of irreversibility from a 
Hamiltonian point of view. 

• The dissipation theory is now formulated with a Hamiltonian: we can quantize it. 
We discussed the concept of quantum decoherence and derived quantum master 
equations for the reduced density matrix. The master equations used here are valid 
in the regimen of weak system-bath coupling, compared with the self-correlation 
times of the bath. 

• We showed that the master equations are thermodynamically consistent, that is, that 
they have as stationary solution the reduced distribution obtained with quantum 
statistical mechanics. We used the opportunity to discuss the differences between 
equilibrium in the classical and quantum cases. In particular, in contrast to the clas- 
sical distribution, the quantum reduced distribution depends on the damping. We 
concluded the equilibrium properties discussing how to properly define and compute 
thermodynamical quantities. 

• With the phase-space formalism we closed the circle. We transformed the equations 
for the density matrix g, into equations for the Wigner distribution. Working in 
phase space and taking the limit h ^ 0, the quantum equations reduce to their 
classical counterparts, explicitly relating the classical and quantum theories of dis- 
sipation. 

• We discussed the continued-fraction technique to solve master equations (obtainment 
of stationary solutions and the linear and nonlinear responses). In spite of the 
demanding pre-treatment and specificity, when it does work it happens to be a 
quite useful method. In particular, it very suited to solve the Caldeira-Leggett 
master equation in phase space for nonlinear potentials. It also allows to study spin 
problems with S" ^ 1, a limit very demanding with other techniques. 
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• We applied the formalisms and the methods of the first part of the thesis to the 
problem of a damped particle in a periodic potential. The main results are: 

— We obtain the Wigner stationary solutions, and we see how the wave nature of 
quantum mechanics takes reflection in the negative zones in the "distribution" , 
persisting at long times. 

— Computing the currents in a cosine potential, we see how increasing the damp- 
ing wc get closer and closer to the classical limit. 

— We tackled the problem of directional motion in potentials without space- 
reflection symmetry. We computed the first quantum corrections to the rectified 
current. The deviations from the classical limit can be understood appealing 
to two complementary quantum phenomena: tunnel and wave reflection. 

— We obtained a estimation of the decoherence times and the inter-well jumps 
from the calculation of the linear susceptibility. 

• Finally, wc addressed the equilibrium and dynamical properties of quantum super- 
paramagnets. The main results are: 

— We obtained dissipative corrections to the Brillouin magnetization law. In 
particular, we could address the entanglement between system and bath at 
very low temperatures from the equilibrium magnetization. 

— We fully characterized the relaxation mechanisms in these systems: overbarrier 
dynamics, and intrawell processes in the potential wells. This allowed us to 
obtain a simple formula for the longitudinal linear response and closed form 
expressions for the relaxation times for any value of the spin S. 

— We generalized our formalism to the nonlinear response. We showed how to 
detect tunnel by means of the nonlinear susceptibility. 

— We produced a simple expression for the transverse susceptibility, which ac- 
counts for the absorption spectra in the ranges of parameters relevant for ordi- 
nary experiments. 
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